EXPERIMENTAL DESIGN AND DATA ANALYTICAL METHODS FOR DETECTING 
AND CHARACTERIZING INTERACTIONS AND INTERACTION THRESHOLDS ON 
FIXED RATIO RAYS OF POLYCHEMICAL MIXTURES AND SUBSETS THEREOF 



This invention was made using funds from grants from the Environmental Protection 
Agency having cooperative agreement number CR827208-01-0 and the National Institute of 
Environmental Health Sciences grant number T32 ES07334-01A1. The United States 
government may have certain rights in this invention. 

DESCRIPTION 

BACKGROUND OF THE INVENTION 

Field of the Invention 
The invention generally relates to the analysis of interactions among agents in a 
combination of agents. In particular, the invention provides statistical methods for determining 
how many (if any) of the agents in the mixture interact, whether the interaction is synergistic or 
antagonistic, and the threshold concentration of agents at which interactions are observed. 

Related Art 

Potentially noxious chemical agents (e.g. metals, manufacturing byproducts, pesticides, 
etc.) are widely distributed in the environment, both as a result of naturally occurring processes 
and as a "side effect" of advancing technology. Humans are exposed to such agents on a daily 
basis, and the protection of human health from the adverse effects of such exposure is crucial. 
Thus, effects of these agents on human health have been and continues to be the focus of intense 
study. 

In the past, studies have typically focused on the effect of a single agent. However, it is 
now recognized that exposure to a single agent is the exception rather than the rule. 
Environmental sources of potentially harmful agents typically contain several or even hundreds 
of chemical substances, so that a person encounters many compounds at once. For example, 
mixtures of several pesticides are generally applied to agricultural cropland, and typical 
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hazardous waste sites may contain mixtures of several hundred different compounds. Thus, it is 
of prime importance to detect, characterize, and predict the impact of mixtures of chemical 
agents on human health. 

Likewise, the availability of a plethora of drugs for treating health-related disorders, 
5 while providing substantial health benefits, often also results in the prescription of multiple drugs 
for an individual. This is often the case with older individuals who tend to have multiple health- 
related issues. Unfortunately, the side effects of interactions among multiple drugs when taken 
together are frequently unknown and unpredictable. 

In the risk assessment of mixtures of chemicals or drugs, a common default assumption is 

10 that at low doses/concentrations, outcomes that are observed are the result of the simple additive 
effects of the individual components of the mixture. However, it is also recognized that at higher 
concentrations, the toxic side effects of chemical agents may be modified by concurrent exposure 
to other chemical agents. 

The detection of chemical interactions, and in particular, the identification of 

1 5 concentrations or conditions under which additivity is not observed, is of special interest. The 
classical method for detecting and characterizing departures from additivity between 
combinations of drugs or chemicals is the isobologram. The isobologram, introduced as a 
graphical tool, is a plot of a contour of constant response of the dose-response surface associated 
with the combination superimposed on a plot of the same contour under the assumption of 

20 additivity. For a two-component mixture, the analysis of an isobologram compares the observed 
isobol (e.g., combination ED 5 o) to the line of additivity. The line of additivity is formed by 
joining the ED50 associated with each of the individual components calculated from the dose 
response data for the individual components. If the isobol is below the line of additivity, a 
synergism is claimed. On the other hand, if the isobol is above the line of additivity, an 

25 antagonism is claimed. However, there are shortcomings associated with the use of 

isobolograms. For instance, the method used in the construction of an isobologram typically 
does not take data variability into account. Additionally, since it is a graphical method, 
isobolograms effectively are limited to the study of combinations of two or three drugs or 
chemicals because of the practical limits involved with data generation and representation. 
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The interaction index, introduced by Berenbaum, provides a convenient method to 
determine and characterize departures from additivity for a combination of c > 2 or 3 
components. The interaction index, n, is defined by 



5 II ~h 2 _|_ _|_ c (1} 

where c is the number of components, Ai, ^2,- • -^c are the doses in combination 
associated with a desired effect, and ED^^CHEMj), i-l,. . .,c is the dose of the z'th component 
that, when administered alone, produces the same effect. When the interaction index n, defined 
in equation (1), is equal to 1 the c components interact additively; when II is greater than 1 the 

10 components interact antagonistically; and when II is less than 1 the components interact 
synergistically. Again, it should be noted that the individual component dose-response 
information is required to calculate the interaction index. The interaction index is directly 
related to the isobologram, i.e., when n=l, the isobol is coincident with the line of additivity; 
when II > 1, the isobol bows above the line of additivity; and when II < 1, the isobol bows below 

1 5 the line of additivity. An advantage of using the interaction index over the isobologram is that 
the interaction index is not limited to combinations/mixtures of just two or three components. 
However, the biological variability associated with the data is not taken into account by the 
interaction index. 

Statisticians frequently use models of the form 

20 

i=i i=i ;=i «=i >i a=i 

t<J i<j<k 



to approximate the relationship between a mean response of interest, ju, and concentrations of c 
chemicals (xi^C2,..-^c) where g(\x) is a user-specified link function. 
25 Carter et al. (reference 4 of Example 6) have shown that a relationship exists between the 

interaction index proposed by Berenbaum and the parameter in a statistical model that is 
associated with the interaction of the components of the combination. Without loss of generality, 
consider that the combination/mixture of interest involves two chemicals and that the response is 
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continuous. Therefore, following the logic of Carter et al., for the linear models case (which is 
readily extended to the generalized linear model case), the relationship between the response and 
the doses or concentrations of the components in combination can be expressed as 

V = Ao + A*i + A 2 *2 + A 2 *12> (2) 

where 

fj. = is the mean response, E(Y), 
/?o is the unknown intercept, 

/?i is the unknown slope parameter associated with the first component, 

is the unknown slope parameter associated with the second component, 
Pn is the unknown parameter associated with the interaction of the two components, 
and x\ and xi are the doses of the respective chemicals. 

From the model defined in equation (2), the EDm/CHEMf) for the respective 
components can be derived to be 

M ~ A) 



ED mfl {CHEM{) : 



A 



ED mfJ (CHEM 2 ) = £- 



02 

15 Thus, after algebraic manipulation, the model defined in equation (2) becomes 

A*i + $2*2 + #2*12 = j 
M ~ Ao V ~ Ao M - Po 

or 

1 P\2 X \2 *\ *2 

~ T^Ao ~ £ " A>)/ + - A))/ ' 
/A /A 2 

Therefore, it follows that when J3\z=0, the combination of components 1 and 2 is additive, i.e., 
20 the isobologram is coincident with the line of additivity and the interaction index equals 1 . 
Similarly, when J3u>0, the combination of components 1 and 2 is synergistic, i.e., the 
isobologram bows below the line of additivity and the interaction index is less than 1; and when 
/?i2<0, an antagonism is present, i.e., the isobologram bows above the line of additivity and the 
interaction index is greater than 1 . This demonstrates the algebraic equivalence between the 
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statistical model and the interaction index. Gennings et al (reference 16 of Example 6) 
demonstrated the experimental convergence of the statistical modeling approach and the 
interaction index. The number of components that can be considered in the statistical model can 
be generalized to c, and data variability is appropriately accounted for in the resulting inference. 
5 The various methods used to determine and classify departures from additivity described 

above utilize both single-compound data as well as combination data. However, such 
methodology does not cover the situation in which single-compound dose-response data are not 
available. For example, as the number of chemicals in a mixture of interest increases, it becomes 
less likely that an investigator will be able to perform an experiment large enough to support the 
10 estimation of these parameters as well as those associated with the dose-response relationship so 
that the potential advantages cannot be realized. 

Up until the present invention, the prior art has failed to provide adequate methodology to 
address these problems. 

1 5 SUMMARY OF THE INVENTION 

The present invention provides methods for analyzing the interaction of the agents in a 
combination of a large number of agents. The analysis is carried out with respect to a 
quantifiable outcome of interest that occurs as result of exposure to the combination of agents. 

20 For example, the present invention provides methods that can determine whether synergistic, 
antagonistic, or no interactions occur among the chemicals in a mixture of chemicals with 
respect to an outcome of interest (e.g. a health disorder) that is caused by exposure to the 
mixture. Alternatively, the methods may be used to analyze interactions of drugs in a mixture of 
drugs that are administered to a patient. The methods of the invention can be used to determine: 

25 whether and how many of the agents in the mixture interact; in certain applications, which of the 
agents are interacting; and the threshold concentration at which the agents begin to interact. In 
one embodiment of the invention, the analysis methods utilize data obtained with individual 
agents of the mixture (e.g. "single chemical data") in addition to the mixture data. However, in 
another embodiment, the analysis methods are carried out in the absence of data generated with 

30 single agents, i.e. it uses only mixture data. Both methods adequately take into account data 
variability. 
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The invention is especially useful for analyzing combinations of a large number of agents 
for which classical methods of analysis would require a prohibitively large number of data 
points. For example, using the methods of the present invention, it is possible to detect 
interactions in a mixture of dozens or even hundreds of chemicals with an economically feasible 
5 study. 

It is an object of this invention to provide a method of detecting an interaction among 
agents in a group using a fixed- ratio ray design, and to determine whether subsets of the agents 
also interact. The method comprises the steps of a. determining an additivity model from single 
chemical data; b. fitting a mixture model in terms of total dose to mixture data from fixed-ratio 

10 rays; c. statistically comparing the additivity model to the mixture model, wherein a difference 
between the additivity model and the mixture model indicates an interaction among agents in the 
group; d. removing a subset of agents from the group; e. repeating steps b and c for agents 
remaining in the group after removal of the subset; and f. determining whether or not the 
remaining agents interact with the subset of agents by utilizing statistical methods based on 

1 5 algebraic manipulations relating full and reduced ray mixture models. The method may also 

include the step of carrying out steps b and c for the subset of agents. The additivity model may 
be graphically represented as an additivity curve and the mixture model may be graphically 
represented as a mixture curve in terms of total dose. Further, simultaneous confidence bands 
maybe determined on the difference between the additivity and mixture curves, or between 

20 mixture curves on full and reduced rays. 

The invention further provides a method of detecting, in a group of agents, the number of 
agents that interact, and determining whether subsets of the agents also interact. The method 
includes the steps of a. fitting a suitable polynomial to experimental data obtained with a 
combination of said agents; b. statistically identifying higher order terms of the polynomial that 

25 are not equal to zero, wherein the number of agents that interact in said group of agents is equal 
to the degree of the higher order terms that are not equal to zero; c. removing a subset of agents 
from the group; d. repeating steps a and b for agents remaining in the group after removal of the 
subset; and e. determining whether or not the remaining agents interact with the subset of agents 
by utilizing statistical methods based on algebraic manipulations relating full and reduced ray 

30 mixture models. The method may also include the step of carrying out steps a and b for the 
subset of agents. Further, single chemical data may also be utilized. The polynomial may be 
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embedded in a generalized linear model or in a general non-linear model. The method may 
include a step of generating a graphical representation of the polynomial. 

The invention also provides a method of determining an interaction threshold for agents 
in a group. The method comprises the step of generating a generalized linear model or general 
5 nonlinear model that permits estimation of the boundaries between a region of additivity of the 
agents and a region of interaction of the agents, in which the boundaries define the interaction 
threshold. In one embodiment of the method, the region of additivity and the region of 
interaction are determined by the steps of 

a. determining an additivity model from single chemical data; 
10 b. fitting a mixture model that incorporates an interaction threshold parameter in terms of 

total dose to mixture data from fixed-ratio rays; and 

c. statistically comparing the additivity model to the mixture model, in which a region of 
difference between the additivity model and the mixture model indicates a region of interaction 
among the agents, and a region of coincidence between the additivity model and the mixture 
1 5 model indicates a region of additivity among the agents. The method may further include the 
steps of d. removing a subset of agents from the group; e. repeating steps b and c for agents 
remaining in the group after removal of the subset; and f. determining whether or not the 
remaining agents interact with the subset of agents by utilizing statistical methods based on 
algebraic manipulations relating full and reduced ray mixture models. The method may further 
20 comprise the step of carrying out steps b and c for the subset of agents. 

Alternatively, in another embodiment of the invention, the region of additivity and the 
region of interaction are determined by the steps of 

a. fitting single chemical data to an additivity model; 

b. fitting mixture data in terms of total dose to a mixture model that incorporates an 
25 interaction threshold parameter, in which the region of additivity is conditioned on results 

obtained in step a; and 

c. statistically comparing the additivity model to the mixture model, wherein a region of 
difference between the additivity and mixture models indicates a region of interaction among the 
agents, and a region of coincidence between the additivity model and the mixture model 

30 indicates a region of additivity among the agents. The method may further comprise the steps of 
d. removing a subset of agents from the group; e. repeating steps b and c for agents remaining in 
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said group after removal of the subset; and 

f. determining whether or not the remaining agents interact with the subset of agents by utilizing 
statistical methods based on algebraic manipulations relating full and reduced ray mixture 
models. The method may also include the step of carrying out steps b and c for the subset of 
5 agents. 

The region of additivity and the region of interaction may be determined by the steps of 
a. fitting an interaction threshold model parameterized with a polynomial function for regions of 
interaction to experimental data obtained with a combination of said agents, and, b. statistically 
testing whether the interaction threshold parameter is different from zero and identifying higher 
10 order terms of the polynomial that are not equal to zero. The method may further comprise the 
steps of 

c. removing a subset of agents from the group; 

d. repeating steps a and b for agents remaining in the group after removal of the subset; 

and 

15 e. determining whether or not the remaining agents interact with the subset of agents by 

utilizing statistical methods based on algebraic manipulations relating full and reduced ray 
mixture models. Single chemical data may also be utilized in the method. 

The present invention also provides a method of designing experiments that achieve a 
target power associated with a test of additivity. The method comprises the steps of 

20 a. specifying the target power, a significance level, a number of mixture dose groups, and 

a magnitude of departure from additivity in terms of total dose; 

b. setting a candidate total sample size; 

c. formulating a design by expressing a design optimality criterion as a function of the 
target power, the significance level, the number of mixture dose groups, the magnitude of 

25 departure from additivity in terms of total dose, and the candidate total sample size, wherein 

optimal locations of total dose groups and optimal allocations of subjects to the total dose groups 
are determined using a direct search algorithm; 

d. calculating a calculated power associated with the design, and 

e. comparing the calculated power to the target power, 

30 f. repeating steps a-e until the step of comparing shows that the calculated power is equal 

to the target power. The method may also include the step of increasing the total sample size 
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used during the step of repeating if the step of comparing shows that the calculated power is less 
than the target power. Conversely, the method may comprise the step of decreasing the total 
sample size used during the step of repeating if the step of comparing shows that the calculated 
power is greater than the target power. 
5 The present invention also provides software programs for implementing all the 

above-described methods. The software programs cause a computer to carry out the steps of each 
method. 

BRIEF DESCRIPTION OF THE DRAWINGS 
1 0 Figure 1. Observed (asterisk) and predicted threshold additivity responses (solid lines) for (a) 
acephate (b) diazinon (c) chlorpyrifos (d) dimethoate (e) malathion. The threshold additivity 
model is given in (2) and model parameters are given in Table 5 of Example 1. 
Figure 2. Observed (asterisk) and predicted (solid line) mean responses along (a) the full five- 
pesticide fixed-ratio ray (0.040: 0.002: 0.031: 0.102: 0.825) and (b) the reduced fixed-ratio ray 
1 5 where malathion was removed from the mixture and the remaining pesticides are at the same 
relative ratios as given in the full ray (0.2286: 0.0114: 0.1767: 0.5833). The threshold 
additivity model (dotted line) fit using single chemical data is provided for reference. Parameter 
estimates are provided in Table 5 of Example 1. Under additivity, the point estimate for the dose 

threshold is estimates as = 34.69 mg/kg along the full five pesticide fixed-ratio ray 

£Mu«r> 
;=i 

20 and as = 6.07 mg/kg along the reduced fixed-ratio ray. (c) The adjusted fitted 

i u i(reduced) 

dose-response curve for the full ray (dashed) is superimposed with the fitted curve for the 
reduced ray (solid). The threshold additivity model (dotted line) is provided for reference. 
Figure 3. The difference between the fitted models for the mixture data along the ray and that 
predicted under the hypothesis of additivity using the single chemical data and the 95% 
25 simultaneous confidence band on the difference in the concentration effect curves for (a) the full 
five pesticide fixed-ratio ray and (b) the reduced ray where malathion was removed from the 
mixture and the remaining four pesticides are at the same relative ratios as given in the full ray. 
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Figure 4. Observed (asterisk) and predicted (solid line) responses along (a) the full five- 
pesticide fixed-ratio ray (0.040: 0.002: 0.031: 0.102: 0.825) and (b) the reduced ray (0.2286: 
0.01 14: 0.1767: 0.5833) for the generalized linear model, given in equation (4) in Example 1, fit 
using only mixture data. The additivity (or first order) model (dotted line) is provided for 
5 reference. Parameter estimates are provided in Table 3 of Example 1. 

Figure 5. Predicted generalized linear interaction model (dotted line) and additivity (first-order) 
model (solid line) fit using single chemical data and mixture data along the full five-pesticide 
fixed-ratio ray (0.040: 0.002: 0.031 : 0.102: 0.825) where the parameter estimates are provided 
in Table 1 of Example 2. The asterisks represent the sample mean motor activity responses. 

10 Figure 6. Generalized linear model along the reduced ray (0.2268: 0.01 14: 0.1767: 0.5833) 
under the null hypothesis (solid line), under the assumption of additivity (dotted line), and under 
the alternative hypothesis (dashed line) that (a) malathion is involved in second-order 
interactions, (b) malathion is involved in second and third-order interactions, and (c) malathion is 
involved in second, third, and fourth-order interactions.- Assumed model parameters are 

1 5 provided in Table 2 of Example 2. 

Figure 7: Observed and predicted responses from the additivity model with estimates given in 
Table 1 of Example 3 for (A) ACE, (B) DIA, (C) CPF, (D) MAL, and (E) DIM. 
Figure 8: Observed (*) and predicted responses based on the SAR mixture model (solid line) 
and assuming additivity (dotted line) as a function of total dose (mg/kg) for the (A) full ray with 

20 all five chemicals and (B) for the reduced ray with MAL omitted. Figure C includes the 

predicted response under additivity (dotted line) and based on the reduced ray mixture model 
(solid line) without MAL as given in B. In addition, Figure C includes an adjusted curve (dashed 
line) from the full ray mixture data (including MAL). 

Figure 9: Observed (*) and predicted responses based on the SANR mixture model (solid line) 
25 and assuming additivity (dotted line) as a function of total dose (mg/kg) for the (A) full ray with 
all five chemicals and (B) for the reduced ray with MAL omitted. Figure C includes the 
predicted response under additivity (dotted line) and based on the reduced ray mixture model 
(solid line) without MAL as given in B. In addition, Figure C includes an adjusted curve (dashed 
line) from the full ray mixture data (including MAL). 
30 Figure 10. Predicted generalized linear interaction model (dotted line) fit with the log link 

function using single chemical data and mixture data along the full five-pesticide fixed-ratio ray 
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(0.040: 0.002: 0.031: 0.102: 0.825) and the additivity (first-order) model (solid line) where the 
parameter estimates are provided in Table 1 of Example 4. The asterisks represent the sample 
mean motor activity responses. 

Figure 11. Generalized linear model along the reduced ray (0.2286: 0.01 14: 0.1767: 0.5833) 
5 under the null hypothesis (i.e. no interactions due to malathion) (solid line), under the 

assumption of additivity (dotted line), and under the alternative hypotheses (dashed line) that (a) 
malathion is involved in second-order interactions, (b) malathion is involved in second and third- 
order interactions, and (c) malathion is involved in second, third, and fourth-order interactions. 
Model parameters are provided in Table 2 of Example 4. 

10 Figure 12. Results of 5, 6, 7, and 10 point D s -optimal designs for the hypothesis given in (10), 
where the alternative hypothesis assumes malathion is involved in two way interactions. 
Parameter estimates are provided in Table 2 of Example 4. Generalized linear model along the 
reduced ray (0.2286: 0.01 14: 0.1767: 0.5833) under the null hypothesis (solid line), under the 
assumption of additivity (dotted line), and under the alternative hypothesis (dashed line) are also 

15 provided. 

Figure 13: Observed versus predicted responses for the fit of the single chemical data. 
Figure 14: Observed and predicted values using the NHEK cell line for (A) arsenic, (B) 
chromium, (C) cadmium, and (D) lead. 

Figure 15: Predicted concentration effect curves along total concentration for the fixed ratio 
20 mixture ray for the curve under the hypothesis of additivity (solid line) and the fitted curve of the 
mixture data alone (dashed line). The asterisks represent the observed sample mean responses at 
each of the mixture points. 

Figure 16: The difference between the fitted models for the mixture data and that predicted 
under the hypothesis of additivity using the single chemical data (solid line) and the 95% 
25 simultaneous confidence band on the difference in the concentration effect curves. 

Figure 17: Illustrations of Isobolograms for a Combination of Two Drugs/Chemicals. The 
dashed line is the 'line of additivity'; when the isobol bows below the line of additivity a 
'synergism' is claimed; when the isobol bows above the line of additivity an 'antagonism' is 
demonstrated. 

30 Figure 18: Observed Responses and Fitted Curve Under the Additivity Model for the Fit of the 
Mixture Data Using the NHEK Cells. 
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Figure 19: Observed and Predicted Responses for the Higher Order Polynomial Model for the 
Fit of the Mixture Data Using the NHEK Cells. • indicate the design locations of the total dose 
values selected by the A i -optimal design. 

Figure 20: Observed and Predicted Responses for the Higher Order Polynomial Model for the 
5 Fit of the Mixture Data Using the NHEK Cells. (Enlarged) 

Figure 21 : Observed and model predicted SDH responses using a generalized linear model with 
power link where 1=0.5. 

Figure 22: Observed and predicted responses for the chlorination mixture along total dose for 
(A) SDH and (B) the transformed g(m) scale. 
10 Figure 23. Isobologram for the interaction of the hypnotic effects of chloral hydrate and ethanol 
in mice. 

Figure 24. Examples of interaction threshold boundaries in a mixture of two chemicals. 
Figure 25. The interaction threshold between ethanol and chloral hydrate (Elliptical Boundary) 
with number of animals responding (out of 6) st each design point. 
15 Figure 26. Contours of constant response (Elliptical Boundary). 

Figure 27. Observed responses from mixture data, predicted curves from the SAR interaction 
threshold model, and the predicted curves from the zero interaction model on the full ray (A) 
and the reduced ray (B) 

Figure 28. The estimated dose-response curve on the full ray projected onto the reduced ray 
20 under the hypothesis of no effect of malathion compared to the observed dose-response curve on 
the reduced ray using the SAR interaction threshold model. 

Figure 29. Observed responses from mixture data, and predicted curves from the SANR 
additivity model on the full ray (A) and the SANR interaction model on the reduced ray (B). 

25 DETAILED DESCRIPTION OF THE PREFERRED 

EMBODIMENTS OF THE INVENTION 

The present invention provides statistical methods for analyzing data obtained by 
exposing subjects to several different concentrations of a mixture or combination of agents, and 
30 observing and quantitating a result or outcome of the exposure. The methods of the present 
invention allow the determination of whether or not the agents in the combination interact to 
produce the result or outcome, or whether the outcome is merely due to the simple additive effect 
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of the agents in the combination. If an interaction is detected, the methods allow the 
determination of how many of the agents interact, whether the interaction is synergistic or 
antagonistic, and the threshold (minimal) concentration of the combination at which the 
interaction occurs. In some embodiments of the invention, the particular agents that interact are 
5 identified. In addition, the present invention provides methods for designing experiments to 
obtain appropriate data for analysis. 

At all concentrations of a combination that is tested, the relative ratios of the agents in the 
combination remain constant. The design used in the methods is the ray design, which can be 
used to study mixtures of drugs or chemicals at a fixed mixing ratio ("fixed-ratio ray"), but in 

10 which the total dose (concentration) varies. By "total dose" we mean the total amount of 
combined agents per dose, concentration or exposure unit of the mixture. This approach 
simplifies analysis of a study by reducing the dimensionality to a set of 2-dimensional dose 
response curves. Further, the methods of the present invention may also incorporate data from 
multiple rays, i.e. data at specific concentrations at two or more fixed-ratios of the agents under 

15 consideration. 

An example of suitable subject matter for the methods of the present invention is the 
interaction of chemicals in a mixture of chemicals of interest with respect to eliciting a health 
disorder in a test subject. In order to obtain data for analysis, test subjects are exposed to 
increasing concentrations of a mixture of chemicals of interest and some indicator that is 

20 symptomatic of the health disorder is measured at each concentration. The relative ratios of the 
amount of each chemical remain constant at each concentration of the mixture, and the relative 
ratios of chemicals that are chosen for use in the mixture should (in order to allow for the most 
relevant influence) have some useful basis. For example, with pesticides, the ratios chosen for 
the mixture might be those most commonly encountered in the environment. Another example 

25 would be the analysis of the interaction of drugs in a mixture of drugs. In this case, the ratios 
might be chosen to mimic typical doses that would be administered to a patient in a "cocktail". 
Some very low concentrations at which the effect is not detectable should be included if possible, 
and the concentration should be increased until the effect is clearly present in the test subjects, or 
until a reasonable, useful conclusion can be drawn from the data. Appropriate control subjects 

30 should be included. 

Using the methods of the present invention, the data so obtained can be analyzed in order 
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to determine whether the symptoms/outcomes that are observed are the result of simple additive 
contributions of each agent (e.g. a chemical) in the mixture, or whether some of the chemicals in 
the mixture interact, resulting in a non-additive outcome. If an interaction is detected, the 
methods of the present invention can be used to determine whether or not the interaction is 
5 synergistic (i.e. two or more of the chemicals interact to bring about a stronger response to the 
mixture than would be predicted if a simple additive effect was being observed) or antagonistic 
(i.e. two or more of the chemicals interact to attenuate the response that would be predicted to 
occur in the absence of interactions), and how many of the chemicals interact. By acquiring 
additional data in which one or more of the chemicals in the fixed ray is removed from the 
10 mixture (i.e. by forming a subset and utilizing a reduced ray), it may also be possible to 
determine which of the chemicals in the subset participates in the interaction, or at least to 
narrow the possibilities. It is also possible to assay the effect of subsets of agents on other subsets 
by analyzing each separately and then in combination in a mixture containing all elements of 
each subset. 

1 5 Further, multiple rays may be analyzed by testing different mixtures of the same agents in 

which the relative rations of the agents are varied among the different mixtures. In other words, 
for a given mixture, the ratio of chemicals within the mixture remains constant at each tested 
dose, but the ratio varies from mixture to mixture (i.e. from ray to ray). 

In addition, using the methods of the present invention, the threshold or minimum 

20 concentration at which an interaction among the chemicals is observed can be identified. 

Determination of such a threshold can be very useful in, for example, recommending guidelines 
for maximal safe exposure to toxic chemicals. Or, is such an interaction threshold is discovered 
for a combination of drugs, it would be of benefit to develop methods directed to administering a 
drug in combinations either above or below the threshold, depending on the desired outcome. 

25 The methods of the present invention allow an analysis of potential interactions in these and 
other systems that will occur to those of skill in the art, and the present invention is intended to 
encompass all such systems. 

The invention also provides methods of experimental design that result in the acquisition 
of data that is particularly useful for the present analytical methods, as described below. 

30 A major advantage of the methods of the present invention is that, because a mixture of 

the agents is tested, it is not necessary to acquire data with each agent individually or in all 
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possible combinations. Thus, it is possible to carry out a meaningful analysis of many agents 
(e.g. dozens or even hundreds) without the necessity of carrying out an impossibly large number 
of single agent experiments. The invention is especially applicable to the analysis of 
combinations of three or more agents for which classical methods of analysis would require a 
5 prohibitively large number of data points. 

The present invention provides two general methods for analyzing data obtained as 
described above. One method, termed Single Agent Required (SAR), utilizes data acquired with 
a mixture in conjunction with similar data acquired individually with single components of the 
mix. The second method, the Single Agent Not Required (SANR) method, requires only data 

1 0 acquired with a mixture/combination. 

In order to carry out both the SAR and SANR methods, a suitable system for 
investigation must be identified, and a database must be generated or located. With respect to 
generating or identifying data for analysis, those of skill in the art will recognize that many 
suitable test systems exist for tabulating data points regarding the effect of a mixture on a 

1 5 particular outcome in test subjects. For example, many cell culture and animal model systems 
exist for ascertaining the result of exposure to mixtures of agents. Examples include but are not 
limited to: the assessment of cell viability in culture after exposure to the agents, or the 
assessment of a particular growth characteristic of the cells (e.g. slow or rapid growth, 
production of a certain product such as an enzyme, induction of mRNA synthesis, transformation 

20 rate, etc); for animal models observable results include but are not limited to viability of the 
animal; development or recession of tumors; the development of a trait such as loss of appetite, 
lack of motor control, increase or decrease in activity, etc., as well as changes at the level of 
tissues, organs or cells (e.g. changes in mRNA synthesis, or change in appearance or functioning 
of tissue, deposition of substances such as cholesterol, etc.). In addition, the outcome that is 

25 detected may be less concrete, e.g. an increased or decreased perception of well-being may be 
observed. Those of skill in the art will recognize that experimental data suitable for analysis by 
the methods of the present invention may be purposefully generated for such analyses, or may 
already have been obtained for some other reason, e.g. in a laboratory or clinical study, during 
the compilation of actuarial tables, etc. Further, the "test subjects" need not be cells or 

30 experimental animals, but can be any entity of interest. Examples include but are not limited to 
persons (e.g. those in clinical trials or studies, or any persons for which data is available), plants, 
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insects, and other living organisms that can be exposed to the mixture of agents, and for which a 
measurable outcome can be recorded. 

In the practice of the present invention, the "agents" in the combination of agents that is 
analyzed maybe of many types. For example, in some cases the agents are chemical entities and 
5 the combination is simply a mixture of the chemical entities. Examples include but are not 
limited to pesticides, drugs, food additives and other ingredients in food (either artificial or 
naturally occurring, e.g. vitamins, minerals, carbohydrates, lipids, etc.), chemical residues from 
water treatment, chemical byproducts from manufacturing processes, metals, etc. However, those 
of skill in the art will recognize that in other cases, one or more "agents" of the combination may 

1 0 be something other than a chemical substance. For example, an agent may be a condition, 

characteristic, phenomenon or activity that the individuals under study experience, engage in or 
are exposed to. Examples of such conditions, characteristics, phenomena and activities include 
but are not limited to receiving a dose or doses of radiation, exposure to sunlight, consumption 
(or lack thereof) of particular food groups or other food items, time spent in an activity (e.g. 

15 meditating, engaging in exercise), height, weight, income level, age, etc. Any "agent" that can be 
quantified or for which a "dose" can be described, may be utilized in a combination of agents 
analyzed by the methods of the present invention. For example, methods of the present invention 
may be used to analyze data in which the variable is a traditional dosage such as number of mg 
per day; or a description of an activity such as "swims (bikes, walks, etc.) a particular number of 

20 times per week"; or variations in income level. In sum, the methods of the present invention may 
be used to analyze any data type that relates a response variable to levels of an independent 
variable that can be placed into fixed ratios. Further, some categorical type variables which are 
not readily described using fixed ratios (e.g. gender, race, adult vs. juvenile, and the like), may 
nonetheless be analyzed by the methods of the present invention. This is accomplished by 

25 designing individual rays as described above and obtaining data for each of the variables of 
interest. For example, the effect of several concentrations of a mixture of chemicals could be 
studied in adult individuals and compared to similar date obtained with juvenile individuals. 

Likewise, the outcomes of exposure to the mixtures of agents that are analyzed by the 
present invention can be of a wide variety. For example, the outcome that is measured may be a 

30 classical clinical manifestation, and may be generally perceived as "negative" (i.e. undesirable) 
or "positive" (i.e. desirable) with respect to health. Examples of a negative outcome of exposure 
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to a mixture of agents include but are not limited to development of a disease condition (e.g. 
cancer, high blood pressure, and many others that will occur to those of skill in the art) or 
physical impairments (e.g. loss of motor control). In contrast, an outcome may be considered 
"positive" e.g. reduction in tumor size, increase in life span, weight loss, pain reduction, etc. 
5 Further, the outcome may be generally perceived as "neutral" or "arbitrary", depending on the 
context, e.g. number of hours spent sleeping. The type of data that is analyzed by the methods of 
the present invention maybe of any suitable type, including but not limited to categorical, count, 
continuous, time to response, etc. 

The SAR method of the present invention may be used to analyze data such as that 

10 described above that has been obtained using several different concentrations or "doses" of a 
mixture/combination of agents along a fixed ratio ray. In the case of SAR, it is also necessary to 
obtain corresponding single chemical data (i.e. individual data for each chemical in the mixture). 
In order to carry out the SAR method, the single chemical data are fit to an additivity model, and 
the mixture data are fit to a mixture model along each ray. Fitting of the single chemical data to 

1 5 an additivity model provides a predicted and estimated graphical representation of the data that 
would be observed if increasing doses of all the chemicals were administered as a mix, and if the 
effect of each chemical in the mix is additive and not influenced by the presence of the other 
chemicals. For example, when using a generalized linear model, the mixture data associated with 
a fixed ratio are fit to a mixture model in terms of total dose: 

20 g(MmiJ=0o + Oi't 

where g(ju mix ) is a specified link function, /? 0 is an unknown parameter associated with the overall 
intercept, and Oi is the unknown parameter associated with the slope relating g(fimb) to the total 
dose, t. By "specified link function" we mean a specified monotone and differentiable function 
relating the mean to total dose in a linear model. Examples of specified link functions include but 

25 are not limited to the logit link often useful for binary data, the log link for count data, and 
identity link for continuous data, and complementary log-log link for percentage data, etc. In 
order to achieve adequate fit to the mixture data, higher order terms in total dose are added as 
necessary. Once both fits are completed, the results are compared. A difference between the 
mixture data fit compared to that of the fit of the additivity data hypothetical mixture data based 

30 on the single chemical data is indicative of a departure from additivity, i.e. the chemicals/agents, 
when in combination, do not elicit the same additive result that is predicted to occur based on 
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data obtained with individual chemicals/agents. Details of such an analysis are given in the 
Examples below, where Examples 1-4 describe both SAR and SANR methods, and Example 5 
describes SANR in particular. 

In addition, it may be desirable to carry out further analyses with a "reduced fixed-ratio 
5 ray", a mixture in which one or more of the chemicals in the original mix (i.e. a subset of the 
original mix) is removed from the mix while the relative ratios of the retained components 
remain the same as in the full ray. For example, if, upon elimination of one or more components 
from the mix, additivity is restored, then this is evidence that the omitted chemical(s) is 
associated with departure from additivity of the mixture. Methods are included in the invention 

10 that allow for statistical hypothesis testing regarding the interaction of subset(s) with the 

remaining agents in the mixture. Such determinations are carried out using statistical methods 
based on algebraic manipulations relating full and reduced ray mixture models. 

In contrast, for the SANR method, single chemical data are not necessary for the analysis. 
This is a very valuable technique for cases where single chemical data may be unavailable or 

1 5 unreliable and cannot be used to describe the additive "no interaction" case. The SANR method 
is based on the assumption of a general parametric form of the underlying response surface such 
that polynomial terms of degree two or greater for a model along a fixed-ratio array are 
associated with interactions among the chemicals in a mixture. Rather than use single chemical 
data, such data are "replaced" by the assumption of the parametric form of the underlying dose 

20 response surface. (However, if single chemical data are available, they can be used to support the 
conclusions obtained with an SANR analysis). In the case of SANR, higher degree terms, which 
are interpreted as being associated with interactions, are included in the model. To account for 
general dose-response shapes (e.g. sigmoidal, exponential, etc.) on the ray in terms of total dose 
caused, for instance, by plateaus, etc., it is necessary to embed the polynomial model of 

25 maximum degree c in a general nonlinear model when a generalized linear model is not 
appropriate. This is done to allow for the polynomial terms to account for the effects of the 
interaction and not the general shape (e.g. sigmoidal, exponential, etc.) of the dose-response 
relationship. Such models will be referred to as "suitable polynomial models". For the case of a 
generalized linear model, a suitable polynomial model is a polynomial model of degree at most 

30 c. In considering the interaction for a mixture of c agents, polynomials with at least 2 and up to 
c-degree terms are considered, i.e. with respect to the term of the polynomial, i, i = 2, . . .c. For 
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example, in a mixture of a full ray of five agents, polynomial models with up to a fifth degree 
term are considered. Upon fitting the full ray data to a polynomial, the significance of the i th 
degree term of the polynomial is interpreted as evidence of an i th degree interaction. By 
"significance of the i' th degree term", we mean that the / th degree term is not equal to zero. 
5 However, the significance of the i th degree interaction is not indicative of which components are 
interacting. In the example of a mix of five agents, ABCDE, if the 4 th degree term is found to be 
significant, then four components are deemed to be interacting. However, five 4-way interactions 
are possible: ABCD, ABCE, ABDE, ACDE, and BCDE. In order to determine which 
combination of four components interacts, it is necessary to obtain data with reduced rays for any 

10 of the possible five groups of four agents where the relative ratios between agents remains the 
same as in the full ray. The details of an SANR analysis are given in the Examples below, where 
Examples 1-4 describe both SAR and SANR methods, and Examples 6 and 7 describe SANR in 
particular. Further, those of skill in the art will recognize that, while no single chemical data are 
required for analysis using the SANR method, if such data are available it may be used in 

1 5 conjunction with an SANR analysis. 

If, as a result of an SAR or SANR analysis, the hypothesis of additivity is rejected, it is 
then of interest to determine total dose regions where interactions occur. To that end, a plot of a 
simultaneous confidence region on the difference in the model fit along the fixed-ratio ray and 
the additivity model over total dose may identify such regions. Total dose levels where the 

20 simultaneous confidence band does not include zero are suggestive of locations and 

characterization of interaction. Through use of simultaneous confidence bands, the user is 
assured of at east the nominal level of confidence. The details of such an analysis are given in 
Example 5 below. From this example, it is clear that along a fixed ratio of the components in the 
mixture, there are regions of synergy, regions of antagonism, and regions of additivity. 

25 In addition, the present invention provides a method to estimate and to test for the 

presence of an interaction boundary in the mixture over all possible ratios of the components, i.e. 
a boundary in the response surface for the combination beyond which interactions among the 
agents in a mixture exist. Example 8 develops methodology for estimating the boundary equation 
which illustrates the methodology with a two-component mixture. Once evidence of departure 

30 from additivity is obtained through the methods developed in this invention, the hypothesis of an 
interaction threshold is of interest. If the null hypothesis of no interaction threshold is not 
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rejected, then the interaction may exist throughout the entire dose region. If the null hypothesis is 
rejected, then the interaction threshold boundary must split the dose region into areas of 
additivity and interaction, which can be identified. Details of such a determination are given in 
Example 2 below. 

5 In addition, the present invention provides an iterative method of designing experiments 

that achieve a target power associated with the test of additivity through the determination of a 
total sample size, determining the location of data points (i.e. total dose values) and allocation of 
subjects to these data points. The location of data points and allocation of subjects are achieved 
through use of optimal experimental design methodology. In order to carry out the method, it is 

1 0 necessary to identify the test system and the agents to be tested. Mixtures/combinations of the 
agents should be designed so that the agents are present in the mixtures in a constant ratio that is 
somehow relevant, e.g. a ratio that is likely to occur in the environment. If c agents are in the 
combination, then c + 1 points (at a minimum) should be placed on a ray in order to optimize the 
criterion of interest associated with the test of additivity. This approach can be applied to tests of 

1 5 the effect of a subset of chemicals on the remaining chemicals (see Example 3). Positioning of 
the points (i.e. dose location) is determined, for example, by A 1 -optimal design or Ds-optimal 
design as illustrated in the Examples. Since statistical power is related to total sample size, 
locations and allocations of experimental subjects to design points (i.e. total doses), the general 
strategy for planning a mixture study is iterative. The basic steps of the design method are: 

20 Step 1 . Specify the target power, significance level, number of mixture dose groups, and the 
magnitude of the departure from additivity in terms of total dose; 
Step 2. Set a candidate total sample size; 

Step 3. Express the design optimality criterion as a function of the values specified in steps 1 and 
2, and determine the location of total dose groups and the allocation of subjects to these groups 

25 that optimize the criterion using a direct search algorithm; 

Step 4. Calculate the power associated with the design determined in Step 3. If the power is 
below the target power, increase the candidate total sample size and repeat Step 3. If the power is 
above the target power, decrease the candidate total sample size and repeat Step 3. Repeat Steps 
3 and 4 until the power associated with the design equals the target power. Details of such design 

30 methods are given in Example 4 below. 
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EXAMPLES 

The following examples are intended to illustrate various embodiments of the invention 
but are intended for purposes of illustration only and should not be construed so as to 
limit the invention in any way. 

Example 1. Detecting Inter action(s) and Assessing the Impact of Component Subsets in a 
Chemical Mixture Using Fixed-Ratio Mixture Ray Designs 

1. INTRODUCTION 

In 1996, the U.S. government passed the Food Quality Protection Act (FQPA), which 
directed the EPA to consider cumulative and aggregate exposure in assessing the risk of 
exposure to chemicals. Due to the variety of agricultural, household, pet, and garden uses of 
organophosphorus (OP) pesticides, the most widely used pesticides in the United States (Aspelin, 
1 994), there is high potential for multiple exposures from multiple routes. Cohen (1 984), 
DuBois (1961), McCollister et al. (1959), and others have studied pesticide mixtures. However, 
these studies have primarily focused on binary chemical combinations and few studies have 
addressed interaction among organophosphorus pesticides. Other authors (e.g., Ma et al., 2002) 
suggest some of these pesticides are human carcinogens. Since humans are exposed to these 
pesticides on a daily bases and there is evidence to suggest exposure leads to adverse human 
health effects, it is of interest to assess the risk of exposure to environmentally relevant 
combinations. 

In studying relevant combinations of OP pesticides, we are interested in detecting and 
characterizing departure from additivity, i.e., determining if exposure to combinations of these 
pesticides results in an increase in the toxic effect above what would be expected under 
additivity. Traditionally, definitions of additivity are based on the classical isobologram (e.g., 
Loewe and Muischnek, 1926; Loewe, 1953). One such definition is Berenbaum's interaction 
index (Berenbaum, 1981), which is algebraically equivalent to the isobologram. The interaction 
index, II, for combinations of c chemicals, is given by 

= 1, additivity ' 

:1, synergism > (1) 
>1, antagonism 

30 where Xi represents the dose of the i th component in combination that yields the desired response 
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(ED m ) and EDf is the dose of the i th component alone that yields the desired response. When II 
equals one, the chemicals are said to be additive; otherwise, departure from additivity is claimed. 
Gennings et al. (1997) illustrated that single chemical dose-response data are necessary and 
sufficient to support the estimation of an additivity model which is algebraically equivalent to the 

5 definition of additivity in ( 1 ) when n= 1 . 

The classical design used to test for departure from additivity in a mixture of c chemicals 
is the factorial design which grows exponentially with c. An alternative to the factorial design is 
the ray design, described by Martin (1942), Mantel (1958), Finney (1964), Brunden and Vidmar 
(1989), and others. Ray designs allow researchers to describe the relationship among multiple 

10 compounds at fixed mixing ratios of interest. For a given ray, the mixing ratio for c chemicals is 

defined by [a\\ a-i. :a c ] where ^Ta f = 1 . Defining t as the total dose along the mixing ray, 

the amount of the i th compound in the mixture is a\t As the number of chemicals, c, in the 
mixture gets large, the experimental effort required to estimate the c+1 -dimensional response 
surface becomes infeasible. Using fractional factorial designs requires the assumption that 

15 certain higher-order interactions do not exist and may result in inaccurate representation of 

nonlinear sigmoid-shaped dose-response relationships. Alternative to focusing on the complete 
response surface, when inference is directed to relevant fixed-ratios of the c chemicals, the 
experimental effort is reduced to include c single chemical dose-response curves and dose- 
response curves in terms of total dose for each fixed-ratio ray of interest. 

20 To motivate the use of these fixed-ratio mixture ray designs consider a study which examines 
exposure to environmentally relevant pesticides. The OP pesticides chosen for study, based on 
usage patterns (i.e. pesticides used on same or similar crops) and market share (i.e. highest 
volume pesticides), were acephate, diazinon, chlorpyrifos, dimethoate, and malathion. Motor 
activity, a count of the number of passes made across a central area, was used to assess toxicity 

25 in healthy adult rats. A full five-pesticide fixed-ratio ray, given by (0.040: 0.002: 0.031: 0.102: 
0.825) for (acephate: diazinon: chlorpyrifos: dimethoate: malathion), was chosen based on 
relative dietary exposure estimates of each chemical as projected by the U.S. EPA Dietary 
Exposure Evaluation Model (DEEM). Since preliminary analysis on single chemical data 
indicated that malathion is not dose responsive with respect to motor activity and malathion 

30 accounts for 82.5% of the total mixture, researchers considered a second fixed-ratio ray, given by 
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(0.2286: 0.0114: 0.1767: 0.5833) for (acephate: diazinon: chlorpyrifos: dimethoate), where 
malathion was removed from the mixture and the four active pesticides were at the same relative 
ratios as those considered in the full ray. For example, the ratio of acephate to diazinon along the 
full ray was 0.040: 0.002 which is a 20:1 ratio. Similarly, the ratio of acephate to diazinon along 
5 the reduced ray was 0.2286: 0.01 1 4 which is also a 20: 1 ratio. 

In what follows, the methods proposed by Gennings et al. (2002) and Meadows et al (2002) 
are extended to consider departure from additivity for threshold models along multiple fixed- 
ratio rays simultaneously. Using the full and reduced rays, methodology is developed for 
assessing the impact of subsets of chemicals (e.g., malathion) on the mixture. Simultaneous 

1 0 confidence bands on the difference between the additivity model and the model estimated along 
a fixed-ratio ray, which may elucidate regions where interactions occur, are also developed. This 
general strategy is developed for two methods of analysis. The first approach, described in 
section 2, is an extension of the work by Gennings et al. (2002), which uses the single chemical 
data for estimation of an additivity model. This model is used to predict the dose-response 

1 5 relationship for a mixture of the chemicals along relevant fixed-ratio rays in terms of total dose 
under the hypothesis of additivity. The observed mixture data are then used to estimate a 
'mixture model' along the ray(s) of interest. The test of additivity is equivalent to the statistical 
comparison of these two curves. If the mixture data support a different relationship than 
predicted from the additivity model estimated by the single chemical data, then departure from 

20 additivity is claimed. This first approach is termed the 'single agents required' (SAR) method of 
analysis. A second approach is also developed which extends the work of Meadows et al (2002) 
who demonstrated that only data along the mixture ray are necessary for detecting departure 
from additivity. In this 'single agents not required' (SANR) method, assumptions are made 
about the underlying relationship along the rays and the significance of higher-order terms in the 

25 polynomial model indicates departure from additivity. In section 3, these methods are extended 
to consider departure from additivity along multiple fixed-ratio rays simultaneously in the 
presence or absence of single chemical data. Tests of hypotheses for making inference about 
interactions due to subsets of chemicals are also developed. The methods are illustrated using a 
mixture of five OP pesticides in section 4. 

30 



23 



2. METHODS: Single Agents Required (SAR) Approach 
2.1. Model Estimation 

In evaluating the risk associated with exposure to mixtures of chemicals, it is often of interest 
to detect a threshold (Schwartz, 1995). Assuming a quasi-likelihood framework (e.g., 
5 Wedderburn, 1974), the threshold additivity model (e.g., Gennings et al., 1997), which relates 
the doses of the chemicals under study to the mean through a link function, g(m), can be 
expressed as 



1 0 g(Madd) is the link function (see McCullagh and Nelder, 1 989) of the response of interest, 
x\ is the dose of the i th single chemical, 
Bo is the unknown parameter associated with the intercept, 
fi j is the unknown parameter associated with the slope of the i th chemical, 
5 is the unknown parameter associate with the threshold, and 

15 8' = ^- is the parameter associated with the dose of the threshold for the i th chemical. 

The model in (2) is described for both an increasing additivity relationship and for a decreasing 
additivity relationship. An additivity model should include mixtures of chemicals that are either 
all nondecreasing or nonincreasing and not a combination of the two cases. If all of the dose 
thresholds are estimated outside of the experimental region resulting in an overparameterized 
20 model, then the corresponding generalized linear model 




(2) 



where: 



is used to replace (2). 
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Let K be the number of fixed-ratio rays under consideration and a (k) = [ai (k ), a 2 (k),- • -,a C (k)] 
define the mixing ratio along the k th fixed-ratio ray. Let JC(k) = [*i(k)> *2(k)>—-,*c(k)] define a vector 

of doses at a given mixture point and t=^] x i(k) define total dose along the k th ray. Following 

1=1 

Gennings et al. (2002) and Meadows et al. (2002), when ray designs are considered, total dose is 
the independent variable along the mixing ray and the amount of the i th compound in the mixture 
along the k th ray is given by a^X. As a result, the threshold additivity model along the k th fixed- 
ratio ray for increasing curves becomes 



8{^add{k)) = 



i 



k if£/w«* 
i i=i 

\fio+t.fi t a m t-S if f,P t a t t*S 

(3) 

Po if*, M f<* 



where 0 l(i) = ^ J P i a i{k) and the parameter associated with the threshold dose under additivity 



1 0 along the k th ray is given by 6[ k) = . Similar results hold for decreasing curves where the 



inequalities in (3) switch direction. If the estimated dose threshold is outside of the experimental 
region, then the additivity model along the kth ray is the corresponding generalized linear model 
given by 

g(M add(k) ) = Po+0 H k) t - 

The model, given in (3), is fit to the single chemical data. Now consider the increasing 
threshold model fit to the mixture data along the k th fixed-ratio ray given by, 

if e' m t<a' ik 

u \(k) l ~ u (k) if 0Kt) tiflt (*) 
where: 

g(Mmix(k)) is the specified link function of the response of interest (McCullagh and Nelder, 1989), 
20 t is the total dose along the fixed-ratio ray, 

Po is an unknown parameter associated with the intercept, 



g(M miX (k)) = 



< 4 > 
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0* m is the unknown parameter associated with the slope along the k th fixed-ratio ray, 

a[ k) is the unknown parameter associated with the threshold along the k th fixed-ratio ray, and 

a ik) = * s parameter associated with the threshold dose along the k th ray. 

In the case where 6[ k) ox a" k) are outside of the experimental region, the threshold model, 
5 given by (3) or (4), is an over parameterized model; in particular, if the dose thresholds have 
estimates less than zero, then there is more than one parameter to describe the intercept. Thus, 
the corresponding generalized linear model is fit. In the case that S' (k) is within the experimental 
region and a" k) is outside of the experimental region, the threshold additivity model is fit to 
single chemical data and the corresponding generalized linear model is fit to the mixture data 

1 0 along the k th fixed-ratio ray. 

A reasonable simplifying assumption is to consider a common intercept for the additivity and 
mixture models. If multiple control groups are experimentally evaluated the assumption of a 
common intercept should be verified by comparing the means of the control groups. As the 
derivatives of the threshold model are not continuous, linearization algorithms may not converge 

1 5 to solutions. Thus, a direct search algorithm (e.g., the Nelder-Mead simplex algorithm; Nelder 
and Mead, 1965) may be used to estimate model parameters using the method of maximum 
quasi-likelihood (McCullagh and Nelder, 1989). 

Before testing for departure from additivity, a goodness-of-fit test is performed on the model, 
given in (4), to ensure that it adequately fits the mixture data. In the event of significant lack-of- 

20 fit associated with (4), higher-order terms are added until an appropriate model is determined. 
Graphical techniques can also be used to assess the fit of the model. Although a goodness-of-fit 
test may indicate that there is no significant lack-of-fit associated with the first order model, a 
plot may indicate that a higher-order model provides a more accurate representation of the data. 
Since the goal is to adequately represent the mixture data along the fixed-ratio rays with a given 

25 parametric model before testing for departure from additivity, both goodness-of-fit tests and 
graphical methods are considered. 
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2.2 Hypothesis Testing 

Following the logic of Gennings et al. (2002), if the relationship between the compounds is 
additive, the parameters along the mixture curves will be equivalent to those predicted under 
additivity. Thus, the overall test for departure from additivity (assuming a common intercept and 
higher-order terms are not necessary to improve the fit of the model) is given by, 

0 m =lM(i 



H 0 :b aid y = 



5 



(5) 



where the pxl vector 7 = [^o.A»Av")^^»^i(i)» a (V^(2)' a (V"'t)' a w^ is a vector of 
model parameters and b a dd is an appropriate 2Kxp contrast matrix. In a more general case where 
higher-order terms are required or the models allow for different intercepts, the additional 
parameters are added to 7 and b a dd is adjusted appropriately. Under additivity, higher-order 
terms are zero and the intercept under additivity is equivalent to those along the fixed-ratio rays. 
A Wald-type test for testing the hypothesis given in (5), is given by 



w= (b7)'[bObT'(bT) 
Mt 



(6) 



where zO is the variance-covariance matrix of 7 and b is any contrast matrix (e.g. b add ). Since 
7 is distributed asymptotically normal (McCullagh and Nelder, 1989) with mean 7 and variance 
til , it follows that W is approximately distributed chi-square with M degrees-of- freedom 
(M=2K for the hypothesis given in (5)). The moment estimate for t is expressed as 

1 „o>,-A-) 2 X 2 



(jv-p)fr {N- P y 

where X 2 is the generalized Pearson statistic (McCullagh and Nelder, 1989), which is 
20 asymptotically distributed chi-square with N-p degrees of freedom. In the quasi-likelihood 
framework, McCullagh (1983) defines the large sample variance-covariance matrix for 7 as 
r[/(7)] _1 , where 1(g) is the expected quasi-information matrix. Let N be the number of 



(7) 
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observations. McCullagh (1983) expressed the expected quasi-information matrix as 

/( 7 ) = D'(V(//))- , D (8) 

where 

Vijly) 0 .» 0 

o run) - o 
; ; o 

0 0 ... V(ju N ) 
5 Replacing y with y in (8), ft = [/(7)]" 1 is a consistent estimate for W (McCullagh, 1983). 
Replacing t with f and W with ft in (6), W is approximately distributed F with M numerator 
degrees of freedom and N-p denominator degrees of freedom, where p is the number of 
parameters. (Note that when higher-order terms are added along the fixed-ratio rays to improve 
the fit of the model or when multiple intercepts are considered, these parameters must be 
1 0 accounted for in the numerator and denominator degrees-of-freedom) . 

In the case that the simultaneous test, given in (5), is rejected, two degree-of-freedom tests 
can be used to detect departure from additivity along individual rays using Hochberg's correction 
for multiple testing, (Hochberg, 1988), which is useful in maintaining the overall significance 
level, a. Hochberg's method is a refinement of the Bonferroni adjustment procedure. It is a 
1 5 sequential procedure in which all p-values are ranked from smallest to largest and the r th p-value 

is compared to , where M is the total number of tests. If p-value^ < — , then 

M-r + 1 w M-r + 1 

the r th hypothesis and all hypotheses with subsequently smaller p-values are rejected. 



and V(fi) = 



2.3. Simultaneous Confidence Band 

20 If the hypothesis of additivity is rejected, it is of interest to determine total dose regions where 
interactions occur. A plot of a simultaneous confidence region on the difference in the model fit 
along the fixed-ratio ray and the additivity model over total dose may identify such regions. In 
order to accommodate the construction of the simultaneous confidence band, differences in the 
two models are considered on the 'linearized' g(u) scale. Consider the simultaneous confidence 

25 band for the k th fixed-ratio ray. For an increasing dose-response curve along the k th fixed-ratio 

ray, define p (k) (t) = g(ju addW )- g(M mix{k) ) as the difference between the transformed response 
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model for the curve under additivity and that fit to the mixture data. Conditioning on values of 
the dose threshold parameters along the ray, S*^ and a" k \ , the threshold models can be written 
as two-phase linear models 

g(^,) = A + (EAW-^ K) , 
>«(")) ' 

where I, c .„ „ is an indicator function that takes on a value of one if the total dose value along the 

(<><>(*)) 

k th ray is greater than the estimated threshold dose under additivity along the k th ray and zero 
otherwise. Similarly / (<>a ^ } is an indicator function that takes on a value of one if the total dose 

value along the k th ray is greater than the estimates threshold dose from the model along the ray 
10 and zero otherwise. 

Let jc correspond to the pxl vector of total dose such that p {k) (t) = x* y. Following Carter 
et al. (1986) and using the Cauchy-Schwartz inequality, it can be shown that conditional on the 
dose threshold parameter values, 

[**'(? -r)j 

pr{x' Qx) 



s Tf F - 



15 This leads to the following conservative probability statement 

U\r-r)) 2 



pr(x* Ox") 

which can be manipulated to yield a conservative simultaneous confidence band on x* 7 such 
that the conditional 100(l-a)% simultaneous confidence band on p (k) (t) is 

P[p L (0 < p (k) (0 < p„ (t) for all f e 91] > 1 - a , (9) 

20 where p L (t) = x'f -[pi(x ' Ox')F p ^_ p .,_J n and p v (t) =X f + [pr{x Qx)F u _ p . x _J 12 . 
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2.4. Assessing the Impact of Subsets of Chemicals 

In addition to simultaneous confidence bands, researchers may be interested in studying the 
effect of particular compounds on the remaining components in a mixture. For example, 
consider a five chemical mixture where researchers are interested in studying the effect of the 
5 fifth chemical on the remaining components of the mixture. This requires the collection of data 
along an additional ray where the fifth chemical is removed from the mixture and the remaining 
components are at the same relative ratios as those given in the full mixture (i.e. 

a,(m . = ). The subscript 'full' denotes values associated with the full five chemical 

a j(full) a j(reduced) 

mixture and the subscript 'reduced' denotes values associated with the mixture where the fifth 
10 chemical was removed. 

Recall for a given dose point t full = x, + x 2 + x 3 + x 4 + x 5 and x 5 - a Hfull) t full . Consider the 
following relationship 

deduced =X l +X 2 +X 3 +X^ 

= _ tfi "'y 5 ^ 'r*** =tfull . (1C 

"~ tfull a S(full)t full 0 ~ a 5(full) ) 

= (1~ a 5(full)) t full 

Given the relationship described above in (10), we can develop hypotheses that compare model 
1 5 parameters along the full and reduced rays. Under the null hypothesis that the fifth chemical 
does not have an effect on the mixture #(%„//;) = g(fi (reduced)) where 

^ ^Kfull^full < a (full) _^ 
/? 0 + 0 X (full$ full - <X(f u U) if fi^full)* full — a {full) 



g(M {f uii)) = 



( l - a 5(full)) 



( X - a 5(full)) 



and 
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SiM(reduced)) ~ 



Pd + ^{(reduced ) f reduced ' a (reduced ) ^ ®\(reduced ) l reduced — a (reduced ) 



Thus, the hypothesis that the fifth chemical does not have an effect on the mixture is given by 

6\(fuii) _ a * 
r - 



Ho interact 7 = 



(!- a 5( full)) 



l (full) ~ u (reduced) 



(11) 



where 7 = [ , # , /? 2 , & , & , /? 5 , 0* l(filll) , a\ m , 0* Kreduced) , a (reduced) ]' and 
1 



0 0 0 0 0 0 0 



5 b in teract = 

is used to test this hypothesis. 



0-10 



(l- a s (/u //)) 

0000000 0 10 -1 



The Wald-type test, given in (6) 



3. METHODS: Single Agents Not Required (SANR) Approach 

10 The general strategy of the SAR method, as specified in the name, is that single chemical 

dose-response data are required to estimate the additivity surface. Investigators may not have the 
resources to conduct single agent studies when the number of chemicals in the mixture is large 
(as in a complex mixture where thousands of chemicals may be present). Or, when the single 
chemical data and the mixture data are not conducted concurrently, there may be evidence of 

15 shifts in the dose-responsiveness of the chemicals over time. Therefore, it is of interest to 
develop a method for testing for departure from additivity for the case where single chemical 
data are not available (so that the additivity surface is not estimable). In addition, we assume the 
mixture data are taken on fixed-ratio rays of the mixture(s) of chemicals. Since we are generally 
interested in risk assessment where threshold models are useful, we develop the SANR method 

20 using threshold models. 



3.1 Test of Additivity 

Carter et al. (1988) demonstrated that when the model parameters associated with interactions 
(i.e., cross-product terms) in a generalized linear model are equal to zero the interaction index is 
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equal to one, which implies an additive relationship. Thus, the significance of cross-product 
terms is associated with departure from additivity, i.e., interaction. 

In the SANR method of analysis the general form of the response surface (interaction) 
model for a mixture of c chemicals is assumed to be given by 

5 gGO = P 0 +IM + £Z ZZZ + Pm..,^x 2 ...x c 

(=1 (=1 j=\ i=l j=\ /=1 

i<j i<j<l 

for a generalized linear model. Following Schwartz et al (1995), the extension of this model, 
which allows for a threshold, is given by 

f J3 0 , A(x)<5) 
8{Mmk) \p o+A {x)-5, A(x)>S] 

A(x) = 2 P,x, +±± P ij x i x j + £ J] PvXfXjXf + + p m ^x x x 2 ...x c 

i-i »=i j=\ i=i ;=i /=i 

i</ i<j<l 

for increasing relationships; the inequalities are switched for decreasing relationships. 
1 0 When ray designs are considered, total dose is the independent variable along the mixing 

ray and the amount of the i th compound in the mixture along the k th ray is given by a^t Let K 
be the number of fixed-ratio rays under consideration and <tyy = [ai(k), fl2(k),- • define the 
mixing ratio along the k th fixed-ratio ray. Let X(k)=[*i(k), *2(k>-— >*c(k)] define a vector of doses at a 

given mixture point and t=]T x i(k) define total dose along the k th ray. As a result, X( k) = a (k) t and 
1 5 the interaction model along the k th ray can be expressed as (Meadows et al., 2002) 
Po if C' + iX*/^-, 



6 m is the unknown parameter associated with the first-order term, 
20 0* (k) for i = 2....c k is the unknown parameter associated with the i th way interaction, 
t is total dose, 
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Ck is the number of chemicals under study along the k th fixed-ratio ray, and 
a* {k) is the threshold parameter.. 

It follows that second order terms are associated with two-way interactions, third order terms are 
associated with three-way interactions, etc. 
5 The model is generally fit simultaneously to the K fixed-ratio rays assuming a common 
intercept. If multiple vehicle control groups are experimentally evaluated, a comparison of the 
control group means is necessary to verify the common intercept assumption. Notice that when 
single chemical data are not available, the parameters in the underlying response surface are not 
estimable, so that the linear parameters 0* (k) are estimated directly and not under the constraint 

10 that 6* m = jj£ /? f fl l(t) . Parameter estimates are found by maximizing the quasi-likelihood 
1=1 

(McCullagh and Nelder, 1989). When the threshold estimate associated with total dose along 
the k th mixing ray is not within the experimental region, the corresponding generalized linear 
model is used, i.e., 

fOW^A.+'Jtt' + iXo'' (12b) 

i=2 

1 5 Following the definition of additivity provided in Section 2 and the methods described by 
Meadows et al. (2002), the parameters associated with interaction along the K fixed-ratio rays, 

namely ( 0 2 * (1) , 0 3 * (1) 0^ , 0 2 * (2) , 0 3 * (2) 0 C * 2(2) , 0 2 * (JC) , 0 3 * w 0* Ck(k) ), are zero under the 

hypothesis of additivity. Define the pxl vector y = [fi Q ,Ot (l) ,Ol m ,...,0* ei(l) ,....,&l (K) ,0* 2(K) , 
0* k(K) ]' as a vector of model parameters and define b at id as a matrix of zeros and ones such 

20 that 

b addT = [#2(1) ' 03(1) >•-> 0 c ,(l) »--.02(*> K t (K) ] • 

An overall hypothesis of additivity, based on the significance of higher-order terms, is given by 
^:b ad dT = 0 (13) 
which can be tested using the Wald-type test described in (6) with scale parameter estimated as 
25 in (7). 
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3.2 Testing the Impact of Subsets of Chemicals 

The methods described in Section 3.1 do not provide specific information about which of the 
compounds are interacting with one another. For example, consider a four chemical mixture 
study where all of the available data are fit to the interaction model given in (12). Let (fli(f U ii): 
«2(fuii): fl3(fuii): a 4 (fuii)) denote the mixing ratio (subscript 'full' denotes the model fit to the mixture 
where the ratios associated with all the compounds are non zero). The model along the fixed- 
ratio ray is given by, 



if 



9\(fidl) ( + X'Vh///' < a (full) 



(14) 



where: 

6\(full) ~ P\ a \(full) + Pl a 2{full) + Pi a l(full) + /?4 a 4(/«//)> 

^2(fuil) = P\2 a \(full) a 2(full) + P\i a \(full) a 3(Jull) + P\4 a \(full) a 4(full) + P 2i a 2( full) 0 \ full) + 

P24 a 2{full) a 4(full) + PiA a 3(full) a 4(fuli)' 
®l(full) ~ P\2l a \(full) a 2(full) a 3(full) + Pm a \(full) a 2(full) a 4(full) + Pm a \(full) a l{full) a 4(full) + 

P2i4 a 2(full) a i(futt) a 4(full)> an ^ 
Q^full) ~ P\7h<\ a \(full) a 2(full) a ?,(full) a 4(full)- 



Suppose that 9* 4{filll) = 0, but 0 3 * (jW/) * 0, i.e. significant third-order interactions are detected. 
Since 0 3 * (/u//) * 0, at least one three-factor interaction exists. Unfortunately, there is not enough 
information on the single ray to determine which of the compounds are involved in this 

15 interaction. Meadows (2001) suggests performing = 4 additional experiments each with a 

combination of three of the compounds at the same relative ratios as those used in the original 

ray (i.e. a ' (/ " /0 = l{reduced) where 'reduced' refers to a mixture where one of the four chemicals 

a j(full) a j(reduced) 

is removed). While these additional experiments allow one to determine which chemicals are 
involved in a particular interaction, the experimental effort may become infeasible as the number 
20 of chemicals under study increases. 
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An alternative to performing 




additional experiments is to consider how a particular 



compound or subset of compounds interacts with the remaining^' components of the mixture. 
Consider the four chemical mixture example described above. Suppose it is of interest to 
determine if the fourth chemical is involved in a three-factor interaction. This can be 
5 accomplished by performing an additional experiment with the remaining three chemicals using 
the same relative ratios as those considered along the original ray. 

Recall the model given in (14) and consider the model fit to the three chemical mixture 
(parameters denoted with subscript 'reduced') given by, 



where: 



\Po if @\(reduced)t + ^ i @i(reduced) t <(X (reduced) 

I i=2 

Po + ^{(reduced)* + J] ® {(reduced)* ~ a (reduced) > ^ ^(reduced)* + J] ^i(reduced)* ^ ° '(reduced) 



(15) 



^(reduced) ~ P\ a \(reduced) + P i a 2(reduced) + Pi a 3(reduced)> 

Ql(reduced) ~ P\Z a i(reduced) a 2(reduced) + P\3 a \(reduced) a 3(reduced) + P 23° 2(reduced) a 3(reduced) > 
@l(reduced) ~ P\23 a \{reduced) a 2(reduced) a 3(reduced)- 

Since the three chemical experiment is performed at the same relative ratios as those in the 
original ray, ^ = ^=±=2- . Therefore, 

a j(full) a j(reduced) 

R a \(fall) a 2(full) a 3(full) _ n a Kreduced) a 2(reduced) a 3(reduced) _ ^(reduced) 
"l23 , .3 -A23 . .3 ~ . .3 

\ a \(full) ) \ a \(reduced) ) ^ a i(reduced) > 

and ^ 3(M) 3 from (14) becomes 

( U l(full)) 
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M 2(/«//) u 3(/^)"4(/«//) 
( a Hfull)) 



91. 



\{full). 

a i(full) a 2(full) a 4(full) 

( a H/ull)f 



^(reduced)) \ u \(full). 
R a 2(full) a 3(full) a 4(full) 

+ P234 : 7i • 

\ a \(full)) 

Now consider the hypothesis given by, 



a i(full) a 3(full) a A{full) 



= 0 



= 0. 



(16) 



L^234j 

In the case that there is not sufficient evidence to reject this hypothesis and the study is 
reasonably powered, we conclude that chemical four is not involved in a three-factor interaction 
and the other three additional experiments suggested by Meadows (2001) are not necessary. If 
there is sufficient evidence to reject the hypothesis, we conclude the fourth chemical is involved 
in at least one three-factor interaction. However, without further experimentation we cannot 
determine which of the three interactions involving this chemical are significant. 

As an alternative to considering the effect of chemical four on three-way interactions alone, it 
may be of interest to test the hypothesis that the fourth chemical is not involved in any 
interactions with the remaining components of the mixture. This is done by simultaneously 
testing the effect of chemical four on two and three-way interactions. The hypothesis that 
chemical four is not involved in any interactions with the remaining three chemicals is given by 



( fl K&«) ) 2 



"urn 



( a \(full)) ( a \(reduced)) 
Where 7 = [fi 0 , 9\ {fill l) , Q 2 (fall) > ^(fall) ' ^reduced) ' &2(reduced) > ^(reduced) ] 



= 0 



(17) 



36 



0 0 - 

0 0 0 



^(reduced) ) 



U l(reduced) ' 

0 0 



1 



In general, consider the c-chemical mixture model given by, 

sO"mix) = A> +o;t+e' 2 t 2 +e\t* +ey + +e]j j * + +9j c . 

First, test H 0 : 6 ■ = 0 Vi = 2, ,c. If there is sufficient evidence to reject this hypothesis, 

additional experiments with j (>h) chemicals, at the same relative ratios as the original 
experiment, will provide further information about the h-factor interaction. The hypothesis of 
interest is given by, 



#0 



(18) 



_( a l(reduced)) H i a i(reduced) )* 

Hq : All h-chemical interactions involving 
1 0 chemicals j+1 , j+2, . . . ,or c are equal to zero. 

(Note that the denominators in the hypotheses given above do not have to be (a {{filll) ) h and 
(^{reduced))'' ■ We could have chosen (a 2{full) ) h and(a 2(reduced) ) h , or (a 3(/ „ U) )* and(a 3(re<W) )* , etc). 
Table 1 provides a list of possible hypotheses when a subset of c chemicals is considered. The 
Wald-type test, given in (6) and (7), is used to test the hypotheses, developed in this section, that 
1 5 consider interactions due to subsets of chemicals. 
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ANALYSIS OF A MIXTURE OF FIVE PESTICIDES 

Summary statistics for the single chemical and mixture data are given in the Tables 2 
and 2a. Single chemical data were collected approximately 18 months prior to the data 
along the foil fixed-ratio ray. The data along the reduced fixed-ratio ray were collected 
5 approximately 6 months following the full ray study. Multiple 'positive control' values 
for the single chemicals were experimentally evaluated in order to verify that the dose- 
response relationship did not significantly shift between the time that the single chemical 
data were generated and the time the mixture data were generated. Preliminary data 
checks were performed to successfully verify that the data across the single chemical and 
1 0 mixture studies were compatible. 



Table 2. Summary statistics for motor activity response associated with single chemical 
experiments. 



Chemical 


Dose 


Mean Motor 


Standard 


Sample Size 




mg/kg 


Activity 
Response 


Deviation 




Acephate 


0 
3 


217.88 
200.13 


35.77 
34.13 


8 
8 




10 


165.88 


25.02 


8 




30 


108.75 


62.51 


8 




60 


58.25 


24.23 


8 




120 


33.25 


27.41 


8 


Diazinon 


0 


206.69 


34.77 


16 




5 


190.88 


28.49 


16 




25 


215.56 


24.02 


16 




50 


183.13 


24.44 


8 




75 


165.69 


33.05 


16 




125 


152.00 


38.65 


8 




150 


76.50 


35.40 


8 




250 


61.25 


47.07 


8 


Chlorpyrifos 


0 


198.56 


24.56 


16 




1 


213.13 


31.42 


8 




2 


192.75 


38.61 


8 




5 


213.37 


31.04 


8 




10 


178.31 


34.02 


16 




20 


157.13 


28.31 


8 




25 


162.00 


32.36 


8 




30 


80.13 


34.85 


8 




50 


49.44 


29.32 


16 


Malathion 


0 


195.86 


19.06 


7 




100 


201.50 


28.38 


8 



41 



500 203.75 28.34 8 



Dimethoate 


0 


195.75 


33.12 


8 




5 


188.25 


24.21 


8 




10 


188.63 


53.05 


8 




25 


107.75 


37.23 


8 




50 


103.75 


51.59 


8 




75 


101.50 


59.57 


8 


Table 2a. Summary statistics for motor activity response associated with mixture data. 


Experiment 


Total 


Mean Motor 


Standard 


Sample 




Concentration 


Activity 


Deviation 


Size 




Dose (mg/kg) 


Response 






Full Ray 


0 


199.43 


20.77 


14 




10 


200.92 


27.10 


12 




55 


167.92 


37.55 


12 




100 


117.08 


47.34 


12 




200 


95.17 


34.03 


12 




300 


72.25 


40.93 


12 




450 


60.08 


46.36 


12 


Reduced Ray 


0 


206.13 


38.24 


8 




1.75 


189.58 


30.82 


12 




9.6 


186.92 


26.87 


12 




17.5 


135.67 


40.46 


12 




35 


84.42 


25.98 


12 




52.5 


43.67 


25.27 


12 




78.8 


48.17 


25.80 


12 



5 

4.1 SAR method of analysis 

The threshold additivity model given in (2) and the threshold mixture model given in 
(4) were simultaneously fit to the single chemical data and mixture data along the two 
fixed-ratio rays allowing for a common intercept. The parameter associated with 
10 malathion was not included in the additivity model as it was not significant. Threshold 
dose estimates from each of the single chemicals were within the experimental region and 
the additivity threshold parameter (5) was significant (p=0.005) (Table 3), which 
indicated a threshold exists under additivity. The dose threshold estimates for the two 
fixed-ratio rays were outside of the experimental region, which indicated that the 
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threshold model was not necessary for the mixture data. Thus, the corresponding 
generalized linear model given by 

s(/w>) 0.W ( 19 ) 

was fit to the mixture data along the two fixed-ratio rays. Tests for the adequacy of the 
5 models fit to the data were performed graphically (Figure 1) for the single chemical data 
and by using the scaled deviance as a goodness-of-fit statistic for the mixture models. 
Plots in Figure 1 indicated that the threshold additivity model provided a reasonable fit 
for the single chemical data. Similarly, a goodness-of-fit test using the scaled deviance 
indicated an adequate fit (p-value of 0.502). 

1 0 The resulting parameter estimates and associated p-values for the additivity and 

mixture models are provided in Table 3. Parameters associated with the full five- 
pesticide fixed- ratio ray are denoted with subscript 'full'. Similarly, parameters 
associated with the fixed-ratio ray where malathion was removed from the mixture are 
denoted with subscript 'reduced'. The slopes associated with the four active pesticides 

15 are all negative and significant (p-value, <0.001), indicating that as the dose of these 
pesticides increases, the mean motor activity response decreases. Similarly, as the total 
dose of the mixture increases along the two fixed-ratio rays the mean motor activity 
response also decreases (p-values of <0.001 along both rays). Point estimates and 95% 
confidence intervals for the threshold doses associated with each single chemical ( 5* ) 

20 are provided in Table 5. Chlorpyrifos has the smallest threshold estimate (4.69 mg/kg) 
and diazinon has the largest point estimate (25.0 mg/kg). All four confidence intervals on 
the dose threshold parameters for acephate, diazinon, chlorpyrifos, and dimethoate do not 
include zero, indicating significant dose thresholds for the four chemicals. Under 
additivity, the plane that connects these dose threshold values is termed the 'additivity 

25 threshold surface'. All dose values below this surface are associated with background 
motor activity responses. Dose values above this surface are associated with motor 
activity responses below background. 

30 
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Table 3. Parameter estimates and resulting p-values for the parameters given in (2) and 
the corresponding generalized linear model to (4) for the five-pesticide additivity model 
and the mixture models along the full and reduced fixed-ratio rays. The estimate for the 



Parameter 


Estimate 


Standard 
Error 


P-value 


Threshold 
Estimate 
(mg/kg) 


95% 
Confidence 
Interval 


Po(add) 


5.303 


0.021 


O.001 






pi (acephate) 


-0.020 


0.002 


O.001 


6.30 


(2.14, 10.45) 


P2 (diazinon) 


-0.005 


4.8E-4 


O.001 


25.0 


(9.54, 40.45) 


p 3 (chlorpyrifos) 


-0.026 


0.002 


O.001 


4.69 


(1.72, 7.66) 


p 4 (Dimethoate) 


-0.014 


0.002 


O.001 


8.72 


(3.15, 14.29) 


5 


-0.124 


0.044 


0.005 






fioffull) 


5.260 


0.042 


O.001 






$0(reduced) 


-0.00306 
5.322 


2.7E-4 
0.045 


O.001 
O.001 






^(reduced) 


-0.023 


0.002 


O.001 







5 NOTE: £fl«iow) = ^ = -0.00307 and £Mo-«i> = T -° 018 - Th 

parameter associated with malathion was removed from the model due to lack of 
significance. 



The hypothesis of additivity using the SAR approach, developed in Section 2, makes 
comparisons between the model predicted under additivity, using single chemical data, 
and the model fit to the data along the fixed-ratio rays. For the data considered here, the 
threshold additivity model, given in (2), was fit to the single chemical data and the 
corresponding generalized linear model, given in (19), was fit to the mixture data along 
the two fixed-ratio rays. Thus, the hypothesis of additivity is defined as 



#o:baddT = 



(20) 

Where 7 =[&(««). A> Pa^^0(M)KmV ^reduced) ^reduced)}' ^ 
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'10 0 0 

0 0.040 0.002 0.031 

10 0 0 

0 0.2286 0.0114 0.1767 



0 -1-10 0 0' 

0.102 0 0 -1 0 0 

0 -10 0-10 

0.5833 0 0 0 0 -1 



This hypothesis was tested using the statistic given in (6) compared to an F^n-io with 
scale parameter given by (7). 
5 Figure 2 provides plots of the observed and predicted responses for the mixture 

data and the predicted threshold additivity model along the full and reduced fixed-ratio 
rays respectively. The means for the mixture data at many of the dose groups along the 
two fixed-ratio rays fall below that predicted under additivity. Results of testing the 
hypothesis of additivity are provided in Table 4. With a p-value of <0.001, there was 
10 sufficient evidence to reject the overall test of additivity, given in (20), and conclude 
significant departure from additivity exists on at least one of the two fixed-ratio rays. 
Using Hochberg's correction for multiple testing, significant departure from additivity 
was detecting along both the full and reduced rays (p-values of O.001 for each). 

15 

Table 4. Overall test of additivity and associated two degree-of-freedom tests for 
significance along the individual fixed-ratio rays. 





Statistic 


Degrees-of-freedom 


P-value 


Overall Test of Additivity 


Hypothesis given in (14) 


7.76 


(4, 467) 


O.001 


Test of Additivity along Individual Rays 


Full Ray 


7.47 


(2, 467) 


<0.001 


Change in intercept 


9.21 


(1,467) 


0.002 


Change in slope 


0.001 


(1, 467) 


0.970 


Reduced Ray 


12.48 


(2, 467) 


O.001 


Change in intercept 


3.23 


(1, 467) 


0.073 


Change in slope 


6.22 


(1, 467) 


0.013 
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Test for the Effect of Malathion 

Hypothesis given in (11) 2.70 (2,467) 0.068 

for generalized linear 

model 

Change in intercept 1.03 (1,467) 0.311 
Change in slope 5.13 (1,467) 0.024 

Figure 3 provides plots of the difference between the additivity and mixture curves on 
the transformed response g(u) and the corresponding 95% simultaneous confidence band 
for the full and reduced rays respectively. Since the curves are decreasing, the difference 
5 was taken as the additivity mean minus the mixture mean; thus, a positive value is 
associated with a more extreme or greater than additive response and a negative 
difference is associated with a less than additive response. The simultaneous bands 
included zero for both the full and reduced rays. Although the test for interaction was 
significant in both rays (Table 4), the simultaneous confidence bands were too wide to 

1 0 elucidate the location of the interaction. 

Finally, since malathion is not dose-responsive alone and it accounts for 82.5% of the 
mixture researchers are interested in determining if malathion has an effect on the 
mixture. This can be achieved by considering differences in the models predicted along 
the full and reduced rays. In other words, if the model fit along the full ray is equivalent 

15 to the model fit along the reduced ray, where malathion was removed and the remaining 
four pesticides are at the same relative ratios as given in the full ray, and there is enough 
power then we can conclude malathion does not have an effect on the mixture. 

The hypothesis for testing for the effect of a particular compound on the remaining 
components of the mixture was developed in Section 2 for the threshold model. For the 

20 mixture data considered here a threshold model was not necessary. Recall under the null 
hypothesis that malathion does not have an effect on the mixture g{m(f u u)) = g{m( re du C ed)) 
where 

S(M(full) ) = fio(fall) + full) 1 full ^ 
_ n ®\(full)Keduced 

~ Po{full) a^a 7 

V 1 a S(fitll) ) 

and 
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S^M (reduced)) ~ Po(reduced) + ^(reduced) 1 reduced ■ 

Thus, the hypothesis that malathion does not have an effect on the mixture is given by 



He, :b. nl 



u a(full) ~ H^reduced) 



0-«3(*B)) 



: = 0< r 



(21) 



5 where t=[P mdd) ,0 x ,P 2 ,fc,P 4 ,S, P %m 6 X{M) 

>Po(reduced)>@\ 



0 0 0 0 0 0 1 
0 0 0 0 0 0 0 



0(reduced)^\(reduced) 

-1 0 
0 -1 



]' and 

. The Wald-type test, given in 



(6) and (7) is used to test this hypothesis. The results of testing the hypothesis of no effect 
due to malathion are provided in Table 4. With a p-value of 0.068, there is marginal 
evidence that malathion has an effect on the mixture. The lack of significance of the 
parameter associated with malathion in the additivity model suggests that the effect of 
malathion is not additive (i.e., malathion interacts with the active pesticides). Figure 2c 
provides the adjusted dose-response curve fit to the five-pesticide mixture overlaid with 
the observed and predicted means along the reduced ray. It is evident from this figure 
that the effect of malathion is a higher dose phenomenon as this is where the solid lines 
and dashed lines are separated the most. 



3.2 SANR method of analysis: Mixture Data Alone 

Although not the case here, for illustrative purposes, suppose that significant shifts 
between the single chemical and mixture studies were observed. This indicates that the 

20 mixture studies are not compatible with the single chemical studies; thus, it is not 

appropriate to combine the data for analysis. In this case, the mixture data alone can be 
used to test the hypotheses of interest. Results of fitting the model, given in (12), 
simultaneously to the mixture data across the two rays are given in Table 5. A model 
building approach was used to consider the level of interaction along the rays. The 

25 fourth-order term along the reduced ray was removed and the third, fourth and fifth 
degree terms along the full ray were removed due to lack of significance. Recall the 
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degree of the term in the model is associated with the degree of the interaction term(s) in 
the assumed underlying response surface model. The hypothesis, given by 

■^0 : & add 7 = \p2(fall) ' ^(reduced) ' ^(reduced) ] = ® 

was used to test the hypothesis of additivity given in (13). Since the p value associated 
5 with additivity was O.001 we may reject the hypothesis of additivity and conclude 
departure from additivity exists along at least one of the fixed-ratio rays. Using 
Hochberg's correction for multiple testing, significant departure from additivity was 
detected along the full and reduced rays (p=0.006 and p=0.002, respectively). With a p- 
value of O.001, we conclude that there are significant third-order interactions along the 
10 reduced ray among the four pesticides; with a p value of 0.006, we concluded that there is 
at least one two-way interaction among the five pesticides together. Plots of the 
additivity and interaction model along each of the rays are provided in Figure 4. 

Table 5. Parameter estimates where data along the two fixed-ratio rays are 
1 5 simultaneously fit to the corresponding generalized linear model for (12). 



Parameter 


Estimate 


SE 


p-value 


Po 


5.313 


0.0365 


O.001 




-0.00497 


0.00073 


O.001 


&2(full) 


5.06E-6 


1.84E-6 


0.006 


^{(reduced) 


-0.00553 


0.00881 


0.530 


@2(reduced) 


-0.00097 


0.000354 


0.006 


®i(reduced) 


0.000010 


3.29E-6 


0.002 


T 


11.70 







Table 6. Overall test of additivity and associated two degree-of- freedom tests for 
20 significance along the individual fixed-ratio rays. 

Statistic Degrees-of-freedom P-value 

Overall Test of Additivity 

Hypothesis given in (13) ~ 21.5 " (3,160) "~ <0.001 
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Test of Additivity along Individual Rays 



Full Ray 


4.759 


(1, 160) 


0.006 


Reduced Ray 


13.02 


(2, 160) 


0.002 


Test for the Effect of Malathion 


Hypothesis given in (14) 


23.2 


(3, 160) 


O.001 



5. DISCUSSION 

The analysis method developed here, which includes the use of single chemical data to 
describe the dose-response relationship assumed under the zero interaction (or additivity) 
5 case, is useful when interest is restricted to making inference along specific rays. 

Appropriate statistical comparisons of the dose-response curve under additivity are made 
to the fitted dose-response curve of the mixture in total dose. Plots of the simultaneous 
confidence bands about the difference in the two models on the g((x) scale may identify 
regions where significant departure from additivity exists. The collection of data along 

1 0 subset rays, where particular compounds of interest are removed from the mixture and the 
remaining components are studied at the same relative ratios as in the full ray, are useful 
in characterizing the effect of the chemicals that have been removed. 

The first analytical strategy presented here has been termed the S AR method, as single 
agent data are required to estimate the additivity surface. By contrast, S ANR method 

15 initially described by Meadows et al. (2002) and Casey (2003) requires the assumption of 
the parametric form of the underlying responses surface. In this case, significant higher- 
order terms in the model indicate departure from additivity. The degree of the largest 
significant term along the fixed-ratio ray indicates the order of the interaction. This 
method of analysis can be used with or without the single chemical data. The choice of 

20 which method to use to detect departure from additivity along fixed-ratio rays is based on 
assumptions the analyst is willing to make and the availability of single chemical data. 

The SANR methods can be used in the presence or absence of single chemical data. 
As illustrated in the example, studies need to be appropriately powered (see Casey, 
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2003). Of particular importance is to consider the reduced power (due to a reduction in 
the total sample size) associated with using only mixture data. We recommend the 
collection of single chemical data when feasible. Considering mixture data alone is 
useful for studies of complex mixtures. 
5 When single chemical data are not available, the parameters associated with the slopes 
of specific chemicals (/?,, J3 2 ,...,ft c ) cannot be estimated unless the number of fixed-ratio 
rays is greater than or equal to the number of chemicals under study. However, reduced 
designs, i.e., designs where the number of fixed-ratio rays is less than the number of 
chemicals under study, provide valuable information about interactions. For example, it 

10 may be reasonable to assume chemicals in the same class have the same slope. In 

addition, constraints can be made across the rays so that the parameter associated with 
additivity,#,* , is common across rays. 

One scenario that can be considered for selecting the same number of rays as 
chemicals is to select a relevant mixing ratio and then, in a systematic order, remove 

1 5 single chemicals, one at a time, keeping the remaining chemicals at the same relative 
ratios as given in the full mixture. Such a design allows for the estimation of slope 
parameters and for the detection of higher-order interactions. In addition, this design 
allows one to consider how a particular compound interacts with the remaining 
components by comparing each of the rays to one another using the methods developed 

20 in Section 3. 

Finally, it should be noted that the hypotheses developed here require the fit of 
higher-order polynomial models. The order of the polynomial follows directly from the 
number of chemicals under study (Meadows et al., 2002). Working with such models 
can cause numerical problems, particularly as the order of the polynomial increases. We 

25 are considering a dose range of 0 mg/kg to 450 mg/kg along the full ray and 0 mg/kg to 
80 mg/kg along the reduced ray. Raising such dose values to powers of two or greater 
results in relatively large numbers; thus, parameter estimates associated with higher-order 
terms become smaller as the order of the polynomial term increases. In some situations, 
such small numbers cause numerical problems. One proposed solution to ease these 

30 numerical problems is to consider scaling the dose range. For example, we could divide 
the doses by 10 or 100 and work with either centigrams or decigrams instead of 
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milligrams. Alternatively, for studies of complex mixtures, it may be reasonable to 
specify the level of interaction researchers are interested in detecting. For example, a 
researcher studying a mixture often chemicals may wish to use the SANR method as a 
screening tool to detect up to five-way interactions. Considering the methods presented 
5 here as a screening tool to detect interactions up to a certain level will ease numerical 
calculation problems and aid researchers in defining specific hypotheses or chemicals of 
interest. For example, researchers may be able to determine that some of the chemicals 
under study are not involved in interactions and, thus, can be removed form the mixture. 
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Example 2. Power and Sample Size Calculations for Determining Whether a Subset of 
Chemicals Interacts with the Remaining Components of a Mixture Using Fixed-Ratio 
Designs 

10 1. INTRODUCTION 

The EPA regulates the use of hundreds of individual pesticides based on acute, 
subchronic, and chronic toxicity assessments conducted in healthy animals (EPA 1999). 
Realistically, however, many pesticides are used in combination or in a pattern that 
results in exposure to multiple pesticides. These pesticides may interact with one another 

15 in such a way that results in unexpected adverse human health consequences. However, 
prior to the Food Quality Protection Act (FQPA 1996) regulatory decisions regarding 
exposure to chemical mixtures often assumed additivity, i.e. zero interaction (see FQPA). 
Regulations that were based on toxicity of these pesticides one at a time may not 
adequately protect human health from the adverse effects of cumulative exposure to 

20 multiple pesticides. 

Although previous studies on pesticide interactions have been performed (e.g., Cohen, 
1984; DuBois, 1961; McCollister et al., 1959), they have primarily focused on binary 
chemical combinations. In addition, few studies have addressed interactions among 
environmentally relevant organophosphorus (OP) pesticides, which are the most widely 

25 used pesticides in the U.S. (Aspelin, 1994). Some of these OPs have been listed as 
priority chemicals for study under the Safe Drinking Water Act amendments. 
Agricultural uses include applications on a variety of food and commercial crops; 
household uses include pet, garden, and home applications. Therefore, there is high 
potential for aggregate exposure, or multiple exposures from multiple routes. Of 

30 particular importance is the detection of interactions, with reasonable power, among 
combinations of multiple pesticides. 

Traditionally, chemical mixtures have been studied using response surface 
methodology (RSM), which is often supported by factorial designs. Although factorial 
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designs adequately allow for the detection and characterization of interactions, the 
number of experimental groups and, thus, the number of observations required becomes 
infeasible as the number of compounds under study increases. In an effort to reduce the 
amount of experimental effort associated with the study of chemical mixtures, tests for 
5 departure from additivity using relevant fixed-ratio ray designs have been proposed by 
Gennings et al. (2002), Meadows et al. (2002a), and Casey (2003). In this case we are 
only interested in making inference along the specific rays of interest, as opposed to 
methods which use designs that require more experimental effort to support the 
estimation of a response surface over a larger region of interest. 

10 Fixed-ratio ray designs were chosen to study interactions among combinations of 

multiple pesticides. The pesticides chosen for these studies, based on usage patterns (i.e., 
pesticides used on the same or similar crops) and market share (i.e., highest volume 
pesticides), were acephate, diazinon, chlorpyrifos, dimethoate, and malathion. Motor 
activity, a count of the number of passes made across a central area, was chosen to assess 

1 5 the toxicity of these pesticides in healthy adult rats. The fixed-ratio ray under study was 
chosen based on the relative dietary exposure estimates of each chemical as projected by 
the U.S. EPA Dietary Exposure Evaluation Model (DEEM) and is given by (0.040: 
0.002: 0.031: 0.102: 0.825) for acephate, diazinon, chlorpyrifos, dimethoate, and 
malathion respectively. 

20 Single chemical data indicate that malathion does not appear to be dose responsive 
with respect to motor activity (see Table 1). This lack of activity along the dose-response 
curve and evidence of interactions along the full fixed-ratio ray (described in Section 4) 
provided motivation for examining a subset ray. This second fixed-ratio ray permits the 
study of the effect of malathion on the mixture (Casey, 2003). Although malathion does 

25 not appear to be dose responsive, it is of interest to investigate whether it interacts with 
the active pesticides. Following the logic of Casey (2003), a comparison of the 
interactions observed along the reduced fixed-ratio ray to those observed along the full 
ray permits the detection and characterization of interactions between malathion and the 
remaining four pesticides. 

30 The reduced ray is given by (0.2286: 0.0114: 0.1767: 0.5833) for (acephate, 

diazinon, chlorpyrifos, dimethoate), where malathion was removed from the mixture and 
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the remaining four chemicals are at the same relative proportions as those considered in 
the full fixed-ratio ray. For example, the ratio of acephate to diazinon along the full 
fixed-ratio ray is (0.040: 0.002) which is a (20: 1) ratio. Similarly, the ratio of acephate 
to diazinon along the reduced fixed-ratio ray is (0.2286: 0.01 14) which is also a (20: 1) 
5 ratio. 

The objective of this paper is to describe the power and sample size calculations, 
which are a function of dose location and allocation of subjects per dose group, 
associated with testing the hypothesis that malathion is involved in interactions with the 
four active pesticides. These methods can readily be applied to tests for departure from 
10 additivity as well. The details of the tests for departure from additivity and tests for 

interactions due to subsets of compounds are summarized below in section 2 and further 
details are provided by Gennings et al. (2002), Meadows et al. (2002a), and Casey 
(2003). 

2. TESTING FOR ADDIVITY 

1 5 Meadows et al. (2002a) and Casey (2003) describe an interaction-model-based method 
of analysis where assumptions are made about the form of the underlying response 
surface. In this case, the significance of higher-order terms in the polynomial 
approximation of the dose-response relationship, expressed as a function of the total dose, 
indicates departure from additivity. Meadows et al. (2002a) showed that, with this 

20 approach, only mixture data along a fixed-ratio ray are necessary for detecting departure 
from additivity. Casey (2003) extended these results to detect departure from additivity 
across multiple fixed-ratio rays in the presence or absence of single chemical data. 

In the quasi-likelihood framework, proposed by Wedderburn (1974), where the data 
are assumed to have mean n and variance rV(/i), the generalized linear model relates the 

25 doses of the chemicals under study to the mean through a link function, g(fi). The 
interaction model for c chemicals can be expressed as 

SCO = Po + + Z2>**/*y+ZXIXw^ + + fim...,x 1 x 2 ..JC c . (1) 

M ;=1 7=1 i=l >1 /=1 

i<j i<j<l 

Let K be the number of fixed-ratio rays under consideration for a mixture of c 
30 chemicals and = [ai (k ), a 2 (k), .^(k)] define the mixing ratio along the k th fixed-ratio 
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ray, such that ]T a i(k) = 1 . Let *(k)=[*i(k)> *2(k> • ■ ■■^c(k)] define a vector of doses at a 

given mixture point along the k th ray and t=^ x j(k) define total dose. Given x ( kj = C(k)t, 
the interaction model along the k th fixed-ratio ray can be expressed as 

s(aw>) = A> +o* m t+e^/ +e; {k f + 

5 (2) 
where: 

i<j i<j<i 
g(lUnix(k)) is tne link function specified by the user (see McCullagh and Nelder, 1989), 
&o is the unknown parameter associated with the intercept, 
10 6* (k) is the unknown parameter associated with the first-order term, 

$* (k) for i = 2....c k is the unknown parameter associated with the i th way interaction, 
t is total dose, and 

Ck is the number of chemicals under study along the k th fixed-ratio ray. 
When single chemical data are available, the parameters associated with the slopes of the 
15 individual chemicals are estimable. Thus, the model fit along the k th fixed-ratio ray is 
given by 

g(M*w) = P 0 +P>° m t+-+p c a Ck{k) t+e; (k) t 2 +e; (k /+ • 

(3) 

The K fixed-ratio ray models are generally fit simultaneously assuming a common 
20 intercept. If multiple vehicle control groups are experimentally evaluated, the 

assumption of a common intercept should be verified by comparing the means of the 
control groups. When single chemical data are available, they are used in combination 
with the mixture data to estimate the intercept and first-order coefficients. In the absence 
of single chemical data, all of the model parameters are estimated using mixture data. 
25 Parameter estimates are found by maximizing the quasi-likelihood criterion (McCullagh 
and Nelder, 1989). 
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Following the logic of Carter et al. (1988), which illustrates the relationship between 
the interaction index, proposed by Berenbaum (1981), and the interaction parameters in a 
statistical model, the additivity model, based on the definition of zero interaction, can be 
expressed as 

1=1 

= A, 

(4) 

Based on this additivity model, the parameters associated with interaction along the K 
fixed-ratio rays, namely 

( 02 ( .) > #30) , 0 2 *(2) > ^3(2) 0* (2 ) , , 0\ (K) , 0* 3{K) (*) ) , are zero under the 

hypothesis of additivity. Define the pxl vector 7 = [/? o ,# 1 * 1) ,# 2 * (1 p...,0* (1) ,...., 

0*(jf) ' 02% > » •••> 0c* x (/o ]' as a vector of model parameters and define b a dd as a matrix of 

zeros and ones such that 

b add y = Kd) • 03(1) » • • ■ » 0q(l) » • • • » 02(K) » • • • » 0 Ck ( JC) ] ■ 

An overall test of additivity, based on the significance of higher-order terms, is given by 
(5) 

A Wald-type test for detecting departure from additivity is given by 

w _ (b7)'[b»bT'(bT) 
Mt 

(6) 

where til is the variance-covariance matrix of 7 and b is any matrix of contrasts (e.g. 
b ad d). Since 7 is asymptotically normally distributed (McCullagh and Nelder, 1989) with 
mean 7 and variance tfl , it follows that W is approximately distributed chi-square with 

K 

M = ^ (c ( . -1) degrees of freedom. The moment estimate for t is expressed as 

;'=1 

~_ 1 yiyy-Mi) 2 _ X 2 
\4l 
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(7) 

where X 2 is the generalized Pearson statistic (McCullagh and Nelder, 1989), which is 
asymptotically distributed chi-square with N-p degrees of freedom. In the quasi- 
likelihood framework, McCullagh (1983) defines the large sample variance-covariance 
5 matrix for for 7 as T[I(y)]~ ] , where 1(g) is the expected quasi-information matrix. 
Replacing 7 with y, Q = [I(y)]' x is a consistent estimate for W (McCullagh, 1983). 
Replacing 7 with f and W with Q in (6), W is approximately distributed F with 

K 

M = ^ (c, - 1) numerator degrees of freedom and N-p denominator degrees of freedom, 

where p=M+K+l for the hypothesis given in (5). 
10 In addition to the hypothesis described above, Casey (2003) developed hypotheses for 
detecting interactions among subsets of chemicals. The hypothesis that a chemical or 
subset of chemicals is not involved in interactions with the remaining components of the 
mixture is given by 



#2(/««) 


@2(reduced) 


( a \(full)) 


( a i(reduced) ) 


Own) 


@3(reduced) 


( fl l(«) 3 


( a \(reduced)) 


Om(fult) 


n* 

m(reduced) 


( a Kfull)) m 


( a \(reduced)) 



15 (8) 

where the reduced ray (parameters denoted with subscript 'reduced') represents the 
mixture study where the chemical or subset of chemicals of interest is removed from the 
mixture and the remaining components are studied as the same relative ratios as given in 
the full ray. These hypotheses compare interaction parameters along the reduced rays 
20 (e.g. the no malathion ray) to interaction parameters along the full ray. Parameter 

associated with the full ray (i.e. the ratios associated with each of the chemicals are non 
zero) are denoted with the subscript 'full'. 

For example, consider a mixture study with three chemicals. Let (ai(f U uj: a2(f u iif- 
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ci3(fuii)) denote the mixing ratio. The model along the fixed-ratio ray is given by, 
where: 

^Kfutl) ~ P\ a i(full) + Pz a 2{fuil) + fi3 a 3(full)> 

&2(fiill) ~ fin a \(full) a 2(fall) + fi\3 a \(full) a 3(full) + P 23° 2( fall)** X full)' 
@3(full) ~ Pm a \(fuU) a 2{JUll) a 7,(fuU)- 

Suppose the third-order term is not significant and it is of interest to determine if the third 
chemical is involved in interactions with the remaining components of the mixture. Data 
are collected along a second fixed-ratio ray given by (a induced)'- ^(reduced)) where 

^KMl = ^reduced) mode j fit alQng reduced fay ig givgn 

a j(full) a i(reduced) 

S(M '{reduced)) = Po + ^reduced)* + ®2(reduced) t 



10 where: 



^(reduced) = M(re<W) + Pz^reduced) ^ 
^2(reduced) = Pl2 a i(reduced) a 2(reduced)- 

Since the reduced experiment is performed at the same relative ratios as those in the 
original ray, fi l2 a J^m = p n a J^dfh^ = g^)_ 

( a i(/«H)) ( a \{reduced)) ( a i(reduced)) 

and — 2(/ " /f) becomes 



(«!(/„«))' 



&2{f«ll) _ & 2(reduced) R a \(full) a i(full) R a 2(ful) a 3(full) 

( \ 2 ~ ( \ 2 1 ( \ 2 ^ 23 ( \2 

\ a \(full)t \ a \(reduced)) \ a \{fuil)) \ a i(Jull)) 

02(ful!) ®2(reduced) _ R a i(Jull) a 3(full) „ a 2(fall) a 3(fall) 



, x2 / s2 / x2 " r A'23 , ,2 ' 

\ a i(full)) ^Hreduced)) \ a i(full)) \ a i(Jull)) 

The hypothesis given by, 

&2(full) ^2(re 



= 0, 
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tests for the effect of the third chemical on the mixture in a pair wise manner (i.e., the 
significance of j8)3 and fe). If there is sufficient evidence to reject this hypothesis, we 
conclude that the third chemical is involved in at least one two-way interaction. Further 
experimentation would be necessary to determine specifically which chemicals are 
5 involved in these interactions. If we fail to reject the hypothesis and the study is 
reasonably powered, we can conclude that the third chemical is not involved in 
interactions with the other two chemicals in the mixture. The objective of this paper is to 
provide methodology to power such studies so that it is reasonable to 'accept' the null 
hypothesis if we fail to reject it. 

10 3. POWER AND SAMPLE SIZE ESTIMATION 

Sample size and power calculations are developed for the tests of additivity based on 
parameter estimates. Under the null hypothesis, the test statistic W, given in (6), has an 
approximate central F-distribution. Under the alternative hypothesis, H a :h^ = (b-y)* 
((b7)* * 0 ) where b is an appropriate matrix of contrast (e.g. b aAA ; b in t era ct), W follows an 

1 5 approximate non-central F-distribution with the same degrees of freedom and non- 
centrality parameter, 1, where 

^ = ^((bT)*)'[bnbT 1 ((b7)*). (9) 

The power of the test is given by, 

EI = P[F(M, N - p, X) > F a (M, N -p, X = 0)], (10) 

20 where F a is the critical value of the central F-distribution with M and N-p degrees of 
freedom. Therefore, in order to calculate power, the non-centrality parameter and the 
total sample size must be specified. 

Power and sample size methods condition on specified model parameters necessary to 
calculate the non-centrality parameter. For the interaction-model-based approach 

25 described in Section 2, single chemical data (when available) are used to estimate the 
intercept and first-order coefficients. Similarly, in the case where data along specific 
fixed-ratio rays are available, they are used to estimate higher-order interaction terms 
along their respective rays. Parameter values for the higher-order interaction terms 
associated with fixed-ratio rays where data are not available are specified based on 

30 interactions the researcher defines as biologically meaningful or from previous 
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experience of the investigator. The value of t is specified based on the variation 
researchers expect to see in the data or from the variation that is observed in the available 
data. 

Given t and the vector of expected model parameters, 7 , McCullagh (1983) defines 

the large sample variance-co variance matrix as zfl = T[I(y)]~ l . I(y) is the expected 
quasi-information matrix given by 



/(7) = D'(V(//))-'D (11) 



where 



d/j M.2 N 



V(M)-- 



V(m) 0 ••• 0 
0 V{n 2 ) ••• 0 



: '•. 0 
0 - V(ji N )\ 

10 Recall N is the total sample size and p is the number of parameters. Notice that each row 
of the matrix D, given in (1 1), corresponds to an observation. Therefore, the variance- 
covariance matrix depends on the sample size allocation. The Nelder-Mead Simplex 
algorithm (Nelder and Mead, 1965), in combination with a bi-section algorithm, can be 
used to maximize the power by optimally allocating the mixture observations to the 

1 5 specified dose locations. 

In order to demonstrate the steps used to calculate power and sample size, define 
power as 

n = f(q,jV,fl;r,7,t,a) (12) 

where: 

20 n eX p is the desired level of power given in (10), 
t is the variation parameter, 

ods the significance level (probability of rejecting the null hypothesis when it is true) 
7 is a vector of model parameters, 
t is a vector of total dose locations, 
25 q is a vector of dose allocations such that Zqj=l , 
N is the total sample size, 
W is the expected quasi-information matrix, and 
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Nqs is the number of subjects allocated to the i dose group. 

Thus, power is a function of the model parameters, the variance-covariance matrix, the 
dose locations, the allocation of subjects, and the total sample size. The following steps 
define the process used to make power and sample size determinations: 

1) Pre-specify dose locations (t), model parameters (t and g), desired level of power 
(n exp ), significance level (a), maximum sample size (N max ), and starting values for 
dose allocation (qo). 

2) In the initial run of the algorithm let No = number of dose groups and let Ni as t - 0, 
where Ni as t is the sample size considered in the previous run of the bi-section 
algorithm. 

3) Use the Nelder-Mead Algorithm to determine the optimal allocation that 
maximizes power 

a) Calculate W 0 (7(7) given in equation (11)) 

b) Calculate the non-centrality parameter, 1 0 , given in equation (9) 

c) Calculate n 0 using the expression given in equation (10) 

d) Using the process defined by Nelder and Mead (1965), determine a new 
value of q and repeat steps 3(a)-(c) until q maximizes P for No 

4) Let P max represent the maximum power determined in step 3 

5) Use a bi-section algorithm to increase or decrease No 

a) If P m ax > riexp + s, where e is a suitably small positive number, then 

N = N 0 

b) If P^x <n exp -e then 

N = N 0 

0 2 
N lasl =N 

6) Continue steps (2)-(4) until n max e (n exp - e, n exp + e) or No > Nmax- 



25 
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4. APPLICATION 

The five-pesticide study described in the introduction is used to illustrate the sample 
size and power methods developed in Section 3. Due to the high potential of aggregate 
exposure to OP pesticides, researchers were initially interested in testing the hypothesis 
5 of additivity, given in (5), along the full five-pesticide fixed-ratio ray, given by (0.040: 
0.002: 0.031: 0.102: 0.825) for (acephate: diazinon: chloryprifos: malathon: 
dimethoate), that represents relative dietary exposure estimates. Single chemical data 
were collected and used to estimate an additivity model along the ray (analysis not 
shown). Since, the motor activity response is a count variable it was assumed that 

10 V(m)=m and the log link function (McCullagh and Nelder, 1989) was used to fit the 
generalized linear additivity model described in Section 2. A plot of this additivity 
model, seen in Figure 5, was useful in specifying the active dose range of 0 mg/kg to 450 
mg/kg along the ray. Prior to the sample size and power calculations developed here, an 
equally spaced design with equal allocation and extra dose points in the low dose region 

1 5 was chosen to study interactions along the full ray. Eighty-six animals were considered. 
Fourteen animals were evaluated in the control group. Twelve animals each were 
evaluated at 10, 55, 100, 200, 300, and 450 mg/kg. Summary statistics for the single 
chemical and mixture data are provided in the appendix. 

Since the mixture contains five chemicals, the single chemical and mixture data were 

20 fit to the fifth-order generalized linear model given in (3). Maximum quasi-likelihood 
estimates were found using PROC GENMOD in SAS® (version 8.2) and are provided in 
Table 1 . A plot of the interaction model is provided in Figure 5. The hypothesis given 
by 



@V f„lh 



(13) 

where 7 = [ A» A > A > A > A > A > 0\ m) , 9* XM) , 9* A{m , 0\ {m ]' and 
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b, dd = 



"0 


0 


0 


0 


0 


0 


1 


0 


0 


0" 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 


0 


0 


0 


0 


0 


0 


0 


0 


0 


0 


1 



, was used to test for departure from 



additivity. With a p-value of <0.001, there was sufficient evidence to reject this 
hypothesis and conclude departure from additivity exists along the full fixed-ratio ray. 
The fifth-order interaction term was not significant (p-value, 0.22); thus, it was removed 
from the model. In addition, the parameter associated with malathion, b 5 , was removed 
from the model due to lack of significance (p-value, 0.85) which indicates the lack of 
activity of malathion. The fourth-order interaction term was marginally significant with a 
p-value of 0.07. Although marginal, we may conclude fourth-order interactions exist 
along the full ray; however, without further experimentation, we cannot determine which 
chemicals are involved in these interactions. 



Table 1. Estimated model parameters where the single chemical and mixture data along 



the full fixed-ratio ray are fit to the model given in (3). 



Parameter 


Estimate 


SE 


p-value 


0o 


5.326 


0.020 


O.0001 


j3 1 (Acephate) 


-0.018 


0.001 


O.0001 


/3 2 Diazinon) 


-0.004 


0.0004 


O.0001 


183 (Chlorpyrifos) 


-0.025 


0.002 


O.0001 


04 (Dimethoate) 


-0.0123 


0.001 


O.0001 


0*(2 nd Order Interactions) 


-3E-5 


1E-5 


0.0097 


0* (3 rd Order Interactions) 


1-.5E-7 


7.2E-8 


0.0345 


6l (4 th Order Interactions) 


-1.9E-10 


1E-10 


0.0714 


7 


3.503 







5 

NOTE: YjPi a i(fm~9\vm = -0-0028. The fifth order interaction term was removed 

15 from the model due to lack of significance (p-value, 0.2203). Similarly, the parameter 
associated with the slope of malathion was removed from the model (p-value, 0.8485). 



The existence of fourth-order interactions along the full fixed-ratio ray and the 
lack of activity of malathion provide motivation for studying interactions between 
20 malathion and the four active pesticides. That is, it is of interest to determine if 
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malathion interacts with the other four pesticides even though it is not active alone. This 
can be achieved through studying a reduced ray, given by (0.2286: 0.01 14: 0.1767: 
0.5833) for (acephate: diazinon: chlorpyrifos: dimethoate) where malathion is removed 
from the mixture and the remaining pesticides are at the same relative ratios as given in 
the full ray. Thus, it is of interest to determine the sample size along the reduced ray 
necessary to test, with reasonable power, the hypothesis that malathion does not interact 
with the four active pesticides. This hypothesis is given by 





@2(reduced) 




( a i(reduced) ) 




@i(.reduced) 


K/i,/o) 3 


( a i(reduced)) 




6 \(reduced) 


.(«!(/«//) ) 4 


( a i(reduced) ) 



(14) 



' @3(fall) > ®A(full) > ®2{reduced) ' ^(reduced) ' ^ ^reduced) ) 



0 0 0 0 0 

(0.04) 2 
0 0 0 0 0 0 

0 0 0 0 0 0 



1 

(0.04) 3 
0 



0 
0 

1 

(0.04) 4 



(0.2286) 2 
0 



(0.2286) 3 
0 



(0.2286) 4 



As specified in Section 3, the maximum sample size, significance level, target power, 
model parameter values, dose locations, and starting values for dose allocation must be 
specified in order to calculate power and sample size. Let alpha be 0.05, the target power 

1 5 be 75%, and the maximum sample size be 80 animals. Estimates for the intercept (&), 
parameters associated with the single chemical slopes (/3 1,&2,& 3, /3*), interaction 
parameters along the full ray {0* 2{m ,9l (fi>n) ,dl (full) ), and the variation parameter (t) are 
provided from the analysis performed on the full ray (Table 1). Preliminary mixture data 
along the reduced fixed-ratio ray would be ideal for specifying the additional higher- 

20 order interaction terms ( 9* 2{reduced) , 0* }(reduced) , 0' 4{reduced) )• However, since these data are not 



65 



available the higher-order interaction parameters along the full ray are used as a guide to 
specify similar parameters along the reduced ray under the alternative hypothesis that 
malathion is involved in interactions with the four active pesticides. Since 6* 2(full) , 0* nfiltt) , 
&l {full) , a i (f u ii), and a^educed), have been specified we can define a model along the reduced 
ray under the null hypothesis (i.e., under the assumption that malathion does not interact 
with the remaining pesticides). For example, consider the case where malathion is not 
involved in any two-way interactions. Under this assumption 



2(full) u 2{reduced) _ q 

(. a \(fiill) ) ( a i(reduced) ) 
@2(full) _ @2(,reduced) 
( a \(full) ) ( a \{reduced) ) 



(15) 



K fuin)' 



" * ( a i(reduced) ) ) ~ 



{(reduced) ) ) ~ u 2(reduced) 4 



1 0 Parameter values along the reduced ray under the null hypothesis that malathion does not 
interact with the active pesticides are provided in Table 2. 

Table 2. Model parameters along the reduced fixed-ratio ray under the null hypothesis 
that malathion does not interact with the remaining pesticides and the alternative 
15 hypotheses that (a) malathion is involved in two-way interactions, (b) malathion is 

involved in two and three-way interactions, and (c) malathion is involved in two, three, 
and four- way interactions. 



Parameter 


No 


Two-way 


Two and 


Two, Three 




Malathion 


Only 


Three-way 


and Four- 




Interactions 






way 




5.326 


5.326 


5.326 


5.326 


(3 i (Acephate) 


-0.018 


-0.018 


-0.018 


-0.018 


/3 2 (Diazinon) 


-0.004 


-0.004 


-0.004 


-0.004 


j8 3 (Chlorpyrifos) 


-0.025 


-0.025 


-0.025 


-6.025 


/5 4 (Dimethoate) 


-0.012 


-0.012 


-0.012 


-0.012 


6l (2 nd Order 


-9.8E-4 


-9.5E-4 


-7.84E-4 


-6.664E-4 


Interactions) 










0* (3 rd Order Interactions) 


2.8E-5 


2.8E-5 


2.52E-5 


1.96E-5 


0 4 * (4 th Order Interactions) 


-2.0E-7 


-2.0E-7 


-2.0E-7 


-1.4E-7 



t 3.503 
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NOTE: =C™/) ="0.0157. 

1=1 

Using the no malathion interaction case as a guide, we can specify model parameter 
values under the alternative hypothesis based on cases that are of interest to investigators. 
5 First, we consider the case that malathion is involved in two-way interactions. In this 
case the third and fourth-order interaction terms represent the no malathion interaction 
case (calculated following the process defined in (15)) and the second-order interaction 
term is specified based on changes in the mean values that researchers define as 
biologically meaningful. We also consider the case where malathion is involved in two 

10 and three-way interactions and the case where malathion is involved in two, three, and 
four-way interactions. The interaction model parameter values along the reduced ray for 
these three cases are provided in Table 2. 

In addition to using the parameter estimates under the null hypothesis or no malathion 
interaction case, plots of the interaction models under the alternative hypotheses and 

1 5 tables of total dose values where a given percent change in the mean value occurs 

between the no malathion interaction case and the specified alternative cases are useful in 
defining biologically meaningful interactions. The plots are provided in Figures 6(a-c). 
The additivity model along the reduced ray is also provided to aid in the interpretation of 
the interaction model specified under the alternative hypothesis. Table 3 provides the 

20 doses where a given percent change in the mean response is observed between the no 
malathion interaction case and the models specified under the alternative hypotheses. 
Since the generalized linear model is nonlinear, the differences in the curves are not 
constant. For example, if we are interested in detecting at least a 5% change in the mean 
for the two-way interaction case, the minimum total dose value that is associated with 

25 such a change is 40.33 mg/kg. 

Table 3. Lowest total dose values along the reduced fixed-ratio ray that result in 5%, 
10%, 15%, and 20% changes in mean responses between the model under the null 
hypothesis and the model that assumes (a) malathion is involved in 2-way interactions, 
30 (b) malathion is involved in 2 and 3-way interactions, and (c) malathion is involved in 2, 
3, and 4-way interactions. Assumed parameter values for each of the models are given in 
Table 2. 
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Change in 


Two-way 


Two and Three-Way 


Two, Three, and 


Mean 


Interactions 


Interactions 


Four- Way 








Interactions 


5% 


40.33 


18.37 


15.78 


10% 


56.36 


28.71 


27.04 


15% 


68.25 


42.98 


77.69 


20% 


77.96 


81.88 


81.46 



Thus far, we have analyzed available single chemical data and mixture data along the 
full fixed-ratio ray to develop a hypothesis of interest and to provide some of the 
parameter estimates necessary to calculate power and sample size. We have used the 

5 available information provided by the analysis, in combination with plots, to specify 
additional model parameter values associated with the alternative hypotheses of interest. 
In order to determine the optimal allocation of the 80 available animals to detect 
interactions involving malathion with 75% power at a significance level of 0.05, we need 
to specify dose locations and starting values for dose allocation. Figures 2(a-c) are 

10 helpful in specifying the dose range and dose locations along the reduced fixed-ratio rays. 
The design considered along the reduced ray is similar to the design considered along the 
full ray where the doses are equally space with extra dose groups in the low dose region. 
The total dose values along the reduced fixed-ratio ray are given by 2, 10, 20, 35, 50, and 
80 mg/kg. Eight control observations were pre-specified along the reduced ray. Initially 

1 5 we assume equal allocation across the remaining dose groups. The Nelder-Mead search 
algorithm (Nelder and Mead, 1965) is used to determine the optimal allocation to the six 
total dose groups. 

The results of the power and sample size calculations for the three cases under 
consideration are provided in Table 4. In the case where we assume that malathion is 

20 involved in two and three way interactions, 80 animals are not sufficient to provide the 
desired level of power. In the case where we assume that malathion is only involved in 
two-way interactions, 77 animals and provide 87% power. However, notice that most of 
the animals are allocated to the 50 mg/kg and 80 mg/kg dose groups. These dose groups 
are located where the largest differences between the model under the alternative 

25 hypothesis (i.e., malathion involved in second order interactions) and the model under the 
null hypothesis (i.e., no interactions due to malathion) occur (see Figure 6(a)). 
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Table 4. Power and sample size results for the hypothesis given in (14) where parameter 
estimates along the full fixed-ratio ray are given in Table 1 and assumed parameter values 
along the reduced fixed-ratio ray under the alterative hypotheses are given in Table 2. 



5 





Two-Way 
Interaction 


Two and Three- 
way Interactions 


Two, Three, and 
Four-way 
Interactions 


Power 


87% 


58% 


76% 


Sample Size Allocation 


Control Group 


8 


8 


8 


2 


0 


12 


8 


10 


1 


12 


7 


20 


1 


12 


8 


35 


1 


13 


9 


50 


33 


13 


9 


80 


33 


11 


9 


Total Sample Size 


77 


81 


58 



It is not likely that allocating one or two animals to particular dose groups will be 
appropriate in real-world situations. Thus, we suggest using these results as a guide and 

10 adjusting the sample size and allocation based on experience of the investigator. For 
example, consider the two-way interaction case presented in Table 4. No animals were 
allocated to the 2 mg/kg dose group and only one animal was allocated to thelO, 20, and 
35 mg/kg dose groups. Consider reducing the sample size at 50 and 80 mg/kg and 
increasing the number of observations at the other dose locations. For example, suppose 

1 5 instead of allocating 8, 0, 1 , 1 , 1 , 33, and 33 and animals respectively to the 0, 2, 1 0, 20, 
35, 50, and 80 mg/kg dose groups we consider allocating 8, 6, 6, 6, 6, 22, and 23 animals 
respectively. This modified allocation results in 83% power, as compared to 87% power 
with original allocation, for testing the hypothesis that malathion interacts with the four 
active pesticides. Modifying the original allocation resulted in reduced power; however, 

20 the decrease in power was not significant (i.e., adjusted sample size allocation results still 
provided reasonable power to test the hypothesis of interest). Alternatively, when 
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possible, extra animals can be added to the dose groups with low allocation leaving the 
allocation to the remaining dose groups unchanged. Adding these extra animals will 
result in increased power since the total sample size increased. In either case, before 
running an experiment with adjusted sample size allocations it is important to determine 
5 the power associated with the hypothesis of interest to ensure the new allocations are 
appropriate. 
5. DISCUSSION 

Power and sample size are important issues when designing experiments to test 
particular hypotheses of interest. Since humans are exposed to hundreds of chemicals 

10 though a variety of sources, it is of particular importance to have sufficient power to test 
the null hypothesis of zero interaction. Without sufficient power we may fail to reject the 
hypothesis when it is false. In this case, the claim of additivity is weak and decisions 
based on this weak claim may result in adverse effects. 

Meadows et al. (2002b) discussed sample size and power issues for testing departure 

1 5 from additivity at specific mixture points where the test for departure from additivity is 
based on a comparison of the predicted mean response under additivity to the true mean 
response. The methods proposed here permit sample size and power determination for 
detecting departure from additivity along multiple fixed-ratio rays simultaneously and for 
testing for interactions involving a particular compound or subset of compounds on 

20 interactions in the presence or absence of single chemical data. These parametric based 
power and sample size methods are appropriate for any hypothesis involving linear 
combinations of the model parameters. In order to use these methods the investigator 
must supply the maximum sample size, the target power, the significance level, specified 
total dose locations, and model parameter values under the alternative hypothesis. 

25 The examples determine the sample size allocation necessary to test, with reasonable 
power, the effect of malathion on the four active pesticides. For ease of notation, the 
methods in this paper were developed for the generalized linear model. In evaluating the 
risk associated with exposure to mixtures, it is often of interest to detect a threshold 
(Schwartz et al., 1995). The methods outlined in this paper are readily extended to the 

30 threshold model (Casey, 2003). Observations are heavily allocated to dose groups where 
larger differences occur, as was demonstrated in the generalized linear model example 
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where it is assumed that malathion is involved in two-way interactions. The examples 
also indicate that equal allocation is not always optimal and larger differences between 
the models specified under the null and alternative hypotheses result in increased power. 
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APPENDIX 



Summary statistics for motor activity response associated with the single chemical data. 


Chemical 


Dose 


Mean Motor 


Standard 


Sample Size 




(mg/kg) 


Activity 
Response 


Deviation 




Acephate 


0 


217.88 


35.77 


8 




3 


200.13 


34.13 


8 




10 


165.88 


25.02 


8 




30 


108.75 


62.51 


8 




60 


58.25 


24.23 


8 




120 


33.25 


27.41 


8 


Diazinon 


0 


206.69 


34.77 


16 




5 


190.88 


28.49 


16 




25 


215.56 


24.02 


16 




50 


183.13 


24.44 


8 




75 


165.69 


33.05 


16 




125 


152.00 


38.65 


8 




150 


76.5 


35.4 


8 




250 


61.25 


47.07 


8 


Chlorpyrifos 


0 


190.00 


14.36 


8 



72 





2 


192.75 


38.61 8 




10 


172.13 


17.55 8 




20 


157.13 


28.31 8 




30 


80.13 


34.85 8 




50 


56.38 


29.98 8 


Malathion 


0 


195.86 


19.06 7 




100 


201.5 


28.38 8 




500 


203.75 


28.34 8 


Dimethoate 


0 


195.75 


33.12 8 




5 


188.25 


24.21 8 




10 


188.63 


53.05 8 




25 


107.75 


37.23 8 




50 


103.75 


51.59 8 




75 


101.50 


59.57 8 


Summary statistics for motor activity response associated with mixture data along the full 


fixed-ratio ray. 








Total 


Mean Motor 


Standard 


Sample Size 


Concentration 


Activity 


Deviation 




Dose (mg/kg) 


Response 






0 


199.43 


20.77 


14 


10 


200.92 


27.10 


12 


55 


167.92 


37.55 


12 


100 


117.08 


47.34 


12 


200 


95.17 


34.03 


12 


300 


72.25 


40.93 


12 


450 


60.08 


46.36 


12 
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Example 3. Analysis of Functional Effects of a Mixture of Five Pesticides using a Ray 

Design 

Introduction 

With support from government agencies, industry, and scientific societies 
5 (including the National Institute of Environmental Health Sciences, the American 
Chemistry Council, the U.S. Environmental Protection Agency, the Society for 
Environmental Toxicology and Chemistry, and the Chlorine Chemistry Council) the 
Society of Toxicology recently spearheaded a series of activities intended to advance the 
scientific understanding of environmental mixtures. An expert Working Group was 

10 charged with evaluating the state of the science on environmental mixtures, providing a 
conceptual framework for future mixtures research, and suggesting potential areas for 
empirical and mechanistic experimentation. A resulting white paper was published 
(Teuschler et al, 2002) which outlined three key ideas with extended discussion. The first 
key idea stated that "toxicology experiments on whole mixtures or mixtures components 

1 5 should include doses at or below the no-observed-effect levels [NOAELs/NOELs] for 
individual mixture components. The mixture components that are tested and their 
relative proportions in the mixture also should reflect those seen in environmental 
samples. In addition, the impact of the unidentified materials in the mixtures should be 
considered." 

20 Existing methods for conducting chemical mixture health risk assessments were 

developed to use available experimental animal data as well as the human health effects 
data in the toxicological and epidemiological literature (Teuschler et al, 2002). These 
methods generally rely on default (e.g., dose addition) assumptions whose validity is 
unknown (ATSDR, 2000a; 2000b; 2000c; U.S. EPA, 2000). Recent congressional acts 

25 mandate regulatory evaluations to be based on the toxicities of the mixtures, not just 

toxicities of the components. Recent congressional acts include a focus on mixtures. For 
example, the Food Quality Protection Act of 1996 directs that assessments of pesticide 
safety include consideration of the risk(s) associated with the cumulative effects of 
chemicals that have a common mode of action while the Safe Drinking Water Act 

30 Amendments, also of 1996, requests the development of new approaches for studying 
complex mixtures. Chemical mixtures risk assessment methods fall into the two general 
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categories of whole mixture approaches (in which complex mixtures are evaluated as 
though they are single entities) and component-based approaches (in which the 
interaction of certain individual components in a mixture are considered to estimate 
toxicity of the mixture) (Teuschler et al, 2002). 
5 The toxicity of a mixture depends on the toxicity of the components and how the 

components interact with each other in a dose-dependent way. In some cases the toxicity 
of a mixture may be adjusted by selecting/changing industrial or environmental processes 
that produce the mixture. For example, a chlorination process and an ozonation- 
chloramination process in disinfecting drinking water may result in different levels of 

10 disinfectant by-products in drinking water. Usage and application rates of fungicides, 
insecticides and herbicides depend on many factors. The selection of particular 
pesticides and/or their application limits could be based on knowledge of which 
chemicals combine synergistically. In general, for industrial/regulatory decision-makers 
to select/change these processes, it is important to identify whether or not the components 

15 of the mixture, including the unidentified fraction, are acting additively. 

Practical methods for assessing additivity in mixtures with many components and, 
in particular the impact of the unidentified materials, have only recently been developed 
(Casey et al, 2003a, b) and are not widely available to decision makers. Such methods 
must accommodate economical and practical designs for more than binary and tertiary 

20 mixtures. The methods must be useful for low-dose exposures and reflect the toxicities 
associated with the relative proportions of the mixture as seen in environmental samples 
including the unidentified fraction. An experimental design that satisfies these 
requirements is the 'ray design.' A ray is defined by a fixed mixing ratio of the 
components in a specified mixture. The dose-responsiveness of the mixture is evaluated 

25 in terms of the total dose/concentration of the mixture with the ratio of the components 
fixed. A statistical model that allows for the evaluation of the low-dose region is a 
threshold model (Cox, 1987). If the researcher is willing to assume the possibility of the 
existence of a threshold, then such a model is useful for estimating the threshold dose, 
below which responses are assumed to be not different from background. When used 

30 with a ray design, the analyst evaluates the dose threshold in terms of total dose for the 
exposure relevant mixture. Based on information/data from single components of the 
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mixture, the predicted dose-response relationship under the assumption of additivity can 
be compared to that observed from the mixture. Such a strategy for evaluating mixtures 
has been previously described (Gennings et al, 1997, 2002; Meadows et al, 2002; and 
Casey et al, 2003a,b). 

5 The objective of this manuscript is to demonstrate this strategy with a mixture of 

five pesticides (acephate (ACE), diazinon (DIA), chlorpyrifos (CPF), malathion (MAL) 
and dimethoate (DIM)). The endpoint of interest is neurotoxicity in rats as measured by a 
gait score dichotomized to indicate presence or absence of a gait abnormality (i.e., 
incoordination, loss of balance, etc.). The fixed-ratio ray under study was chosen based 

10 on the relative dietary exposure estimates of each chemical as projected by the U.S. EPA 
Dietary Exposure Evaluation Model (DEEM) and is given by (0.040: 0.002: 0.031 : 0.825: 
0.102) for ACE, DIA, CPF, MAL, DIM, respectively. The data from the five chemical 
mixture study are referred to as 'full ray' data. In order to evaluate the effect of one 
pesticide (MAL) on the dose-response relationship of the other four pesticides, a 'reduced 

1 5 ray' was experimentally evaluated where the remaining four pesticides are fixed at the 
same relative proportions as considered in the full ray, i.e., (0.229: 0.011: 0.177: 0: 
0.583) for (ACE: DIA: CPF: MAL: DM). (Note: Such a strategy can be extended for 
identifying the impact of the unidentified fraction in a complex mixture.) 
Additivity Model 

20 The definition of additivity we use is given by Berenbaum (e.g., 1 985) and is 

based on the classical isobolograms for the combination of two chemicals (e.g., Loewe 
and Muischnek, 1926; Loewe, 1953). That is, in a combination of c chemicals, let £,• 
represent the dose of the i th component alone that yields a fixed response, and let x, 
represent the concentration/dose of the i th component in combination with the c agents 

25 that yields the same response. According to this definition of additivity if the substances 
combine additively, i.e., with zero interaction, then 

C Y 

i=l ^ 

If the left-hand side of (1) is less than 1, then a synergism can be claimed at the 
combination of interest. If the left-hand side of (1) is greater than 1, then an antagonism 
30 can be claimed at the combination. As (1) is the equation of a plane in c dimensions, this 
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definition of additivity implies that under additivity contours of constant response are 
planar. In what follows, the additivity models satisfy both the definition given in (1) and 
the fundamental notion of 'no interaction' (Gennings et al, 2003), i.e., that the rate of 
change (slope) in the response as a function of the i th chemical does not change in the 
5 presence of other chemicals. 

Response surface methodology, often supported by factorial designs, is the 
classical statistical experimental approach for testing for departure from additivity. 
Alternatively, fixed-ratio ray designs have been proposed (e.g., Gennings et al, 2002, 
Meadows et al, 2002) to reduce the amount of experimental effort when the exposure 

10 region of interest is restricted to relevant mixing ratios. Two analysis strategies are 

demonstrated here and result in similar conclusions. The first approach is similar to that 
described by Gennings et al (2002) applied to a threshold additivity model, and is termed 
the 'single agents required' (SAR) method; the second approach (Meadows et al, 2002 
and Casey et al, 2003a) is termed the 'single agents not required' (SANR) method. 

1 5 The general strategy begins by using the single chemical data to fit a threshold 

additivity model. The threshold additivity model is used to estimate the dose-response 
relationship along the fixed-ratio ray(s) of interest (in terms of total dose) under the 
hypothesis of additivity. This additivity model is used to provide power/sample size 
calculations to design the mixture study with adequate power to detect biologically 

20 meaningful interactions. The mixture data are then experimentally generated and fit to a 
threshold model in terms of total dose. For the SAR method, the mixture model is 
statistically compared to the predicted model under additivity. If the models are 
determined to be different, then evidence of departure from additivity is inferred. 
Otherwise, no departure from additivity is inferred. When the studies have the power to 

25 detect a biologically meaningful departure from additivity and there is not significant 
evidence of departure from additivity, additivity can be claimed. 

The SANR method assumes a general parameterization of the underlying 
multidimensional response surface (in this case 6 dimensions). The result is that the 
terms in a polynomial model along a fixed ratio ray are associated with interaction terms 

30 in the underlying model (Meadows et al, 2002; Casey et al, 2003a). For example, a 



77 



significant cubic term in the mixture model is an indication of a three-way interaction 
among the chemicals in the mixture. 

The threshold additivity model (e.g., Gennings et al, 1997) is of the form 



where 

g(/x A ) is the specified link function (McCullagh and Nelder, 1989) of the response of 
interest, 

x\ is the dose of the i th single chemical, 

Po is an unknown parameter associated with the overall intercept, 
Pi is an unknown parameter associated with the slope of the i th pesticide dose response, 
8 is the unknown parameter associated with the threshold for the additivity model. 
Using this model, the parameter associated with the dose of the threshold for the i th 

chemical is given by 6* = — , i=l,. . .,5. 

Pi 

If the model fits a threshold outside the experimental range (resulting in an over- 
parameterized model), then the corresponding generalized linear model is used, i.e., 



Let the mixing ratio of chemicals of interest be denoted as (ai : dq, : a3 : a» : a 5 ) such that 

j£ a, - 1 . Following Gennings et al (2002) and using the definition of additivity given in 
i=i 

(1) and the additivity models in either (2) or (3), the slope along the fixed-ratio ray design 
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under additivity is given by ^ i a j J3 i . 

The biological endpoint of interest is an indicator of neurotoxicity based on evidence of 
the presence or absence of a gait abnormality (Moser, 1995). In order to appropriately 
5 constrain the probability of a response to be between 0 and 1, a logit link function was 
used, i.e., g(/x) = logO /(l- n)). It was assumed throughout that the variance changes as a 
binomial random variable such that Var(Y)= Tfi (1- n) (e.g., McCullagh and Nelder, 
1989). 

Parameter estimates for the threshold model in (2) are found using the maximum 
10 quasi-likelihood criterion (e.g., McCullagh and Nelder, 1989) in a Nelder Mead 
algorithm (Nelder and Mead, 1965). Parameter estimates for the generalized linear 
model given in (3) were found using the maximum quasi-likelihood criterion using a 
Fisher scoring algorithm in Proc GENMOD in SAS (version 8.2). A moment estimate 
for t (McCullagh and Nelder, 1989) was used. Adequacy of the fit of the model to the 
1 5 data was assessed graphically and by comparing the scaled deviance to a xVp distribution 
where N is the total sample size and p is the number of model parameters. 

Initially the threshold model given in (2) was fit to the data. The dose threshold 
estimates were all outside of the experimental region; thus the corresponding generalized 
linear model given in (3) was used as the additivity model. Parameter estimates and 
20 associated p values are provided in Table 1 . Plots of the observed and model predicted 
responses are provided in Figure 7. The single chemical data are adequately represented 
by the additivity model (p=0.341, using scaled deviance as an assessment of fit). 

Table 1 : Parameter estimates and associated p values from the additivity model given in 
25 (3). The slope for MAL was not included as no gait abnormalities were observed in the 



experimental range for MAL. The estimate for the scale parameter was f =1 . 1 87. 



Parameter 


Estimate 


Standard Error 


P value 




-5.130 


0.7061 


O.001 


Pi (ACE) 


0.1868 


0.0368 


O.001 


h (DIA) 


0.0338 


0.0063 


<0.001 


P 3 (CPF) 


0.1605 


0.0268 


<0.001 


MDM) 


0.1711 


0.0373 


O.001 
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In addition to gait score, other biological and neurochemical endpoints were taken. 
Sample size estimates (n=12/mixture group) were calculated for motor activity data to 
detect a 25% decrease (i.e., more negative) in the slope of the dose-response curve along 
5 the fixed-ratio ray from that under additivity with 70% power using a Wald-type test (as 
described in Casey et al, 2003c) with a two-sided test with 5% significance. When the 
decrease in the slope was 30%, the resulting power was 84%. 

As the mixture studies were not run concurrently with the single chemical studies, 
it is important to demonstrate that the single chemical dose-response curves do not shift 

10 across the studies (mixture and single chemical). 'Positive control' values of the single 
chemicals were included in the mixture studies for this purpose. The resulting mixture 
experiment consisted of a vehicle control (n=14), 'positive controls' for ACE alone (6 
and 30 mg/kg; n=6, 14, respectively), CPF alone (10, 20, 30 mg/kg; n=6, 8, 6, 
respectively), DIAZ alone (50 and 125 mg/kg; n=6, 14, respectively), MAL alone (350 

15 mg/kg; n=14), DIM alone (10, 25, 50 mg/kg; n=6, 8, 6, respectively), and six total dose 
mixture concentrations (10, 55, 100, 200, 300 and 450 mg/kg; n=12 rats per dose group 
for a total of 72 rats exposed to a mixture of the chemicals) at a fixed mixing ratio ray of 
(0.040 : 0.002 : 0.031 : 0.825 : 0.102) for (ACE: DIA: CPF: MAL: DIM). Similarly, the 
design for a second fixed-ratio ray was based on a reduced number of the pesticides 

20 (omitting MAL) with the remaining pesticides at the same relative proportions as in the 
'full' mixture study. The experiment consisted of a vehicle control (n=8), positive 
controls (n=8 per group) for ACE (30 mg/kg), CPF (20 mg/kg), DIA (125 mg/kg), and 
DIM (25 mg/kg), and six total dose mixture concentrations (1.75, 9.6, 17.5, 35, 52.5, and 
78.8 mg/kg; n=12 rats per dose group) at a fixed mixing ratio ray of (0.229: 0.01 1 : 

25 0. 1 77: 0.583) for (ACE: DIA: CPF: DIM). Casey et al (2003a,b) describe the motor 
activity data in the mixture of these five pesticides. 
'Single Agents Required' Method 

Following the logic of Gennings et al (2002), for comparison to the additivity 
model along the fixed-ratio rays given by (3), the mixture data were fit to a 'mixture 

30 model' along each ray in terms of total dose. 
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In order to achieve adequate fit to the data, higher order terms in total dose are added to 
(4) when necessary. 

The five vehicle control groups for the single chemical data and the vehicle 
control groups for the two mixture studies resulted in no gait abnormalities. Therefore, 
5 since the background rate is similar across the studies, the additivity model in (3) and the 
mixture model in (4) for both the full and reduced rays were fit simultaneously in an 
overall model with a common background parameter, /3o. In addition, the two mixture 
studies included at least one 'positive control' group for each of the single chemicals 
alone. Preliminary analyses (not shown) of all of the data combined did not find 

10 evidence of a shift in the dose response curves from the original single chemical studies 
using the positive control data (p=0.515). Therefore, the overall model was based on the 
single chemical data and the mixture data combined. 

The overall model included a common intercept term, linear terms for each of the 
single chemicals (excluding MAL), a linear and quadratic term in total dose for the full 

1 5 ray, and a linear term in total dose for the reduced ray (described above), which 

adequately represented the data (p=0.672, i.e., no indication of lack of fit). The resulting 
parameter estimates and associated p values are provided in Table 2. As the linear 
parameters associated with total dose along both the full and reduced rays are significant, 
there is evidence that as the total dose of the mixture increases, the probability of a gait 

20 abnormality also increases. 

Table 2: Parameter estimates and associated p values from the full SAR analysis. The 
estimate for the scale parameter was f =1.44. The overall hypothesis of additivity was 
rejected (p=0.040, with 3 degrees of freedom (df)). The hypothesis of additivity along 
25 the full ray was rejected (p=0.046, with 2 df) but was not rejected along the reduced ray 
(p=0.486, 1 df). The hypothesis of no malathion effect was rejected (p=0.016, 2 df). 



Parameter 


Estimate 


Standard Error 


P value 


Po 


-3.9663 


0.4905 


<0.001 


Pi ACE 


0.1474 


0.0364 


<0.001 


p 2 DIA 


0.0256 


0.0054 


<0.001 


P 3 CPF 


0.1250 


0.0217 


O.001 


p 5 DIM 


0.1335 


0.0328 


<0.001 


Bid) full ray linear 


0.0472 


0.0101 


O.001 
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02(1) full ray quadratic 


-0.0001 


<0.0001 


0.007 


9 K2) reduced ray linear 


0.1166 


0.0209 


O.001 



From the parameter estimates given in Table 2, plots of observed and predicted 
responses along both fixed-ratio mixture rays are provided (Figure 8) based on the 
5 mixture model and under the assumption of additivity. Following the approach in 

Gennings et al (2002), as the curve predicted using the mixture model for the full mixture 
ray falls significantly above the predicted dose response curve under the assumption of 
additivity (Figure 8A) (p=0.046), it can be inferred that the overall effect of the specified 
fixed-ratio of the five pesticides is associated with a greater than additive, or synergistic, 

10 relationship. Further, there is not a significant difference (Figure 8B) in the dose 
response curves from the mixture model and that predicted under the assumption of 
additivity for the mixture (omitting MAL) along the reduced ray (p=0.486). Following 
the logic of Casey et al (2003b), this indicates that the interaction (i.e., synergy) in the 
full ray may be associated with the presence of MAL. Casey et al (2003b) developed a 

1 5 method for combining the results from multiple rays with the same relative ratios (like 
the full and reduced rays here) in the same figure. If the adjusted dose-response curve 
from the full ray is identical to the dose-response curve on the reduced ray, then it 
indicates the chemical or subset of chemicals that were removed in the reduced ray do not 
interact with the remaining chemicals. For these data, testing for the difference in the 

20 corrected and corresponding parameters from the full and reduced rays (Casey et al, 

2003b) indicates that the presence of MAL significantly (p=0.016) changes the shape of 
the dose response curve in the full ray. This is elucidated in Figure 2C where the dashed 
curve corresponds to the adjusted dose-response curve from the five-pesticide mixture 
plotted with the dose response curve without MAL (solid curve). These dose-response 

25 curves are significantly different due to the presence of MAL. The total dose that is 

associated with a 50% response (ED 5 o) for the four pesticides in the absence of malathion 
is about 34 mg/kg. In the presence of malathion it is roughly 19 mg/kg. This nearly two- 
fold difference indicates the magnitude of the effect of the synergism at a 50% response 
level. It is apparent from Figure 8C that the magnitude of the synergism is dependent on 

30 the level of the effect (i.e., 10%, 50%, etc). 
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'Single Agents Not Required' Method 

Following the methods of Meadows et al (2003) and Casey et al (2003a), a 
comparison analysis strategy to the SAR method is based on an assumption of a general 
parameterization of the underlying six-dimensional (five chemicals and one response 
5 dimension) response surface. Meadows et al (2003) showed that polynomial terms of 
degree two or greater for the model along a fixed-ratio ray are associated with 
interactions among the chemicals in a mixture. Casey et al (2003a) generalized this result 
to the case of multiple mixture rays and termed this approach the 'single agents not 
required' (SANR) method. As the name suggests, the method is applicable to the case 

1 0 where single agent data are not available for estimation of the additi vity relationship. 
This information is 'replaced' with the assumption of the parametric form of the 
underlying response surface. However, if single agent data are available, they can be 
used to support the estimation of the corresponding parameters in the model. 

In comparison to the mixture model in (3), here we allow for higher degree terms, 

15 which are interpreted as being associated with interactions. The interaction model for a 
mixture of c chemicals is given by 

iK/0 = /?o+I/V' (5) 

For example, for data from the full ray with five pesticides, a polynomial model with up 
to a fifth degree term was considered; for the data from the reduced ray with four 

20 pesticides in the mixture, a polynomial model with up to a fourth degree term was 

considered. The significance of the i th -degree term is interpreted as evidence of an i th - 
way interaction (i=2, ...,c), but is not suggestive of which i components are interacting. 
For example, in the mixture analysis of five pesticides, (for ease of notation denoted as A, 
B, C, D, E) there are five possible four-way interactions: ABCD, ABCE, ABDE, ACDE, 

25 and BCDE. The significance of the 4 th degree term in the model given in (5) indicates 
that one or more of these four-way interactions are present. Further experimentation is 
needed in order to infer which components are interacting (Meadows, 2002; Casey et al, 
2003a,b). Estimation of model parameters follows similarly to that described in the 
previous section. 

30 The initial model fit to the gait score data had a common intercept term across the 
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studies (5 single chemical and two mixture studies), linear terms for each single chemical 
study for estimating the additivity response, a fifth-degree polynomial for the full ray 
data, and a fourth-degree polynomial for the reduced ray data. All of the data are used to 
estimate the additive part of the model, but only the mixture data were used to estimate 
the corresponding higher-degree terms. From this model, the simultaneous test for the 
significance of the higher-degree terms is equivalent to a test of additivity. For these 
binary gait score data, there was an indication of departure from additivity (p=0.049, with 
7 df). In order to make the model more parsimonious, we used a backward elimination 
approach to delete terms that were not statistically significant. We reduced the model for 
the full ray first while keeping the reduced ray fully parameterized. Once the model for 
the full ray was reduced to having only significant terms, then the model for the reduced 
ray was simplified. None of the higher-order terms for the reduced ray were significant; 
however, we kept the pure quadratic term for flexibility in the model. 

The resulting parameter estimates for the 'final' model, standard errors, and p 
values are provided in Table 3. As all of the linear terms are positive and significant, we 
conclude that as any of the doses of the chemicals increase, the probability of a gait 
abnormality also increases. The overall test of additivity was rejected indicating the 
significance of at least one higher degree term (p=0.041, 3 df). Figure 9 A presents a plot 
of the dose-response curve in terms of total dose for the predicted model compared to 
what it would be under additivity. Starting with the highest degree term, since the fourth 
degree term is significant for the full ray, we conclude that at least one four-way 
interaction exists. Without further information, it is not evident which four chemicals are 
involved. 



25 Table 3: Parameter estimates and associated p values from the final SANR analysis. 
The estimate for the scale parameter was f =1.42. The overall hypothesis of additivity 
was rejected (p=0.021, with 4 df). The hypothesis of additivity along the full ray was 
rejected (p=0.041, with 3 df) but was not rejected along the reduced ray (p=0.175, 1 df). 



Parameter 


Estimate 


Standard Error 


P value 


Po 


-4.244 


0.5486 


<0.001 


pi ACE 


0.1615 


0.0384 


<0.001 


P 2 DIA 


0.0276 


0.0057 


O.001 


P 3 CPF 


0.1343 


0.0235 


<0.001 


PsDIM 


0.1531 


0.0370 


<0.001 
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0 2 (i) fall ray quadratic 


0.000699 


0.000280 


0.013 


63(1) full ray cubic 


-0.00000553 


0.00000245 


0.024 


04(n full ray quartic 


0.0000000105 


.00000000517 


0.042 


0 2 (2) reduced ray quadratic 


-0.0000659 


0.000456 


0.148 



By comparison, none of the higher degree terms are significant along the reduced 
ray (p=0.175). Therefore, there was no evidence of departure from additivity among the 
5 four pesticides (ACE, DIA, CPF, and DIM; see Figure 9B). Recall the relative ratios of 
these pesticides are equivalent to those used in the full ray. Since there is an indication of 
at least one four-way interaction on the full ray where MAL is present, and no indication 
of departure from additivity on the reduced ray where MAL is absent, there is evidence to 
suggest that MAL interacts with at least some or all of the other four pesticides. 

10 Following the work of Casey et al (2003a), the hypothesis of no MAL interaction with 
the other pesticides was rejected (p=0.014, 3 df). Figure 9C depicts the same two curves 
as in Figure 3B with the additional curve (dotted) obtained from the full ray, corrected for 
the proportion of the pesticides that were included in the reduced ray (here 1-.825 since 
MAL was 82.5% of the mixture). The difference in the dotted and solid curves indicates 

1 5 the effect MAL has on the mixture. 
Discussion 

In summary, both analysis strategies resulted in similar conclusions, namely that 
(1) ACE, DIA, CPF and DIM when given alone had a significant effect on a gait 
abnormality and MAL was not dose responsive, (2) there is a significant interaction 

20 among the five pesticides along the fixed-ratio mixture ray which is associated with a 
synergistic effect, (3) there is not a significant departure from additivity among the four 
pesticides (omitting MAL) along the reduced mixture ray, and (4) adjusted comparisons 
across the full and reduced rays indicate that MAL interacts with the other pesticides. 
Casey et al, (2003a,b) describe in more detail the SAR and SANR methods in the 

25 analysis of motor activity in the same rats that were used here for the gait score response. 
Their conclusions for motor activity are similar to those reported here for a gait 
abnormality. 

More generally, the methods described here can be readily generalized to other 
types of chemicals, different numbers of components in a mixture, and to different types 
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of endpoints. For example, the gait score endpoint used here resulted in proportional 
data; motor activity described by Casey et al (2003a,b) is a count endpoint; continuous 
endpoints or time to response endpoints could also be used. The methods based on 
exposure relevant fixed-ratio ray designs allow for dramatic decreases in the required 
5 experimental effort for studying mixtures of many chemicals as compared to the classical 
factorial designs. The design for the study of the mixture of five pesticides essentially 
required adequate support for seven dose-response curves (five single chemicals alone 
and two total dose mixture rays). Positive controls were added in order to demonstrate 
the comparability of the studies not run concurrently. By comparison a full factorial 

10 design for five chemicals with reasonable support of the shape of the dose response curve 
(i.e., not just two dose points per chemical) would have required more experimental 
effort. For example, a full factorial design for a mixture of five chemicals each at three 
levels has 243 dose groups. Of course, fractionated factorial designs are available. But, 
as the number of chemicals in the mixture increases, the experimental effort for even 

1 5 these designs quickly becomes impractical. 

Meadows et al (2002) describe an example where the single chemical data were 
not appropriate for comparison to the mixture data. In such cases, the interaction model 
based methods may be used to estimate departure from additivity with only mixture data 
along a ray to support the estimation of the model parameters. In screening methods for 

20 mixtures of many chemicals, it may not be feasible to supply dose-response curves for 
each component in the mixture. By making the assumption of a parameterization for the 
underlying response surface, analysis of data along a fixed-ratio ray may result in 
hypothesis testing for departure from additivity. 

An important feature to the analysis strategies is the power and sample size 

25 calculations for the mixture study. The single chemical data were initially evaluated 
using an additivity model. Power and sample size methods (Meadows et al, 2003; and 
Casey et al, 2003c) were used to determine dose locations and sample sizes to yield an 
acceptable level of power to detect departure from additivity to a biologically important 
degree. Here a change in the parameter associated with the slope of at least 25% from 

30 that predicted under additivity was considered biologically meaningful. If the study is 

adequately powered and no evidence of departure from additivity is found, then additivity 
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can be claimed. By comparison, if the study is not adequately powered, then lack of 
evidence of departure from additivity may be a power problem and not an indication of 
additivity. 

Using both the SAR and S ANR methods departures from additivity were found 
5 for the given mixture of the five pesticides. This claim of interaction applies to the 

specified mixing ratio and may not be true for other mixing ratios. The characterization 
of the interaction may also change along a fixed-ratio ray. For example, Gennings et al 
(2002) demonstrated in a mixture of four metals that the relationship among the 
chemicals changed from synergistic, to additive, to antagonistic along the specified ray. 

1 0 Casey et al (2003b) developed methods of determining the location of interaction for a 
fixed ratio of chemicals through use of simultaneous confidence regions. 

When studying motor activity with this same combination of pesticides, Casey et 
al (2003a,b) concluded that malathion interacted with the remaining four pesticides. 
Interestingly, the effect did not occur until total doses (of the four pesticides omitting 

1 5 malathion) exceeded about 45 or 50 mg/kg. In the analysis of gait score described here, 
Figures 8C and 9C suggest a malathion effect at much lower total doses. An interaction 
with malathion was not totally unexpected, since a frequently-cited example of OP 
mixtures is the potentiation of malathion toxicity by other OPs, which inhibit the 
carboxylesterase-mediated hydrolysis of malathion (e.g., Murphy and DuBois, 1957). 

20 Given the large amount of carboxylesterases in the body, such a kinetic interaction would 
be expected more at higher dose levels than were observed in the present study. 

The risk assessment process is complicated by the fact that environmental 
exposures frequently may involve mixtures of chemicals rather than a single compound. 
Public concern about such mixtures led to new requirements under the Food Quality 

25 Protection Act (1996) to assess risk of pesticide mixtures that have a common mode of 
toxicity. The methods described here may be useful in providing an experimentally 
feasible way of studying exposure-relevant mixtures instead of regulating chemicals 
based on default assumptions of additivity. In addition, if departure from additivity is 
concluded, the impact of important components of the mixture can be assessed by 

30 comparison of 'full' and judiciously chosen 'reduced' rays of interest. 



87 



References 

ATSDR (Agency for Toxic Substances and Disease Registry) (2000a). Guidance Manual 
for the Assessment of Joint Toxic Action of Chemical Mixtures. 

5 

ATSDR (Agency for Toxic Substances and Disease Registry) (2000b). Interaction 
Profile for: Chlorinated Dibenzo-p-Dioxins, Hexachlorobenzene, p,p'-DDE, and 
Methylmercury. 

10 ATSDR (Agency for Toxic Substances and Disease Registry) (2000c). Guidance for the 
Preparation of an Interaction Profile. 

Berenbaum, M.C. (1985) The expected effect of a combination of agents: The general 
solution Journal Theoretical Biology , 114,413-431. 

15 

Casey, M., Carter, W.H., Jr., Gennings, C, Moser, V., Simmons, J.E. (2003a) A single 
agent not required (SANR) approach for detecting interaction(s) and assessing the 
impact of component subsets in a chemical mixture using a fixed-ratio ray design, 
Journal of Agricultural, Biological and Environmental Statistics , under review. 

20 

Casey, M., Gennings, C, Carter, W.H., Jr., Moser, V., Simmons, J.E. (2003b) A single 
agent required (SAR) approach for detecting interaction(s) and assessing the impact 
of component subsets in a chemical mixture using a fixed-ratio ray design, Journal of 
Agricultural, Biological, and Environmental Statistics , under review. 

25 

Casey, M., Gennings, C, Carter, W.H., Jr., Moser, V., Simmons, J.E. (2003c) Power and 
sample size determination for testing the effect of subsets of compounds on mixtures 
along fixed-ratio rays, Journal of Agricultural, Biological, and Environmental Statistics , 
under review. 

30 

Cox, C. (1987) Threshold dose-response models in toxicology, Biometrics 43: 51 1-523. 

Gennings, C, Schwartz, P., Carter, W.H., Jr., Simmons, J.E. (1997) Detection of 
departures from additivity in mixtures of many chemical with a threshold model, 
35 Journal of Agricultural, Biological, and Environmental Statistics 2: 198-211. 

(2000) Erratum: Detection of departures from additivity in mixtures of many chemical 
with a threshold model, Journal of Agricultural, Biological, and Environmental 
Statistics 5: 257-259 

40 

Gennings, C, Carter, W.H., Jr., Campain, J.A., Bae, D., Yang, R.S.H. (2002) Statistical 
analysis of interactive cytotoxicity in human epidermal heratinocytes following 
exposure to a mixture of four metals, Journal of Agricultural, Biological, and 
Environmental Statistics 7:58-73. 



88 



Gennings, C, Carter, W.H., Jr., Simmons, J.E., Carchman, R., Carney, E. (2003) An 
attempt to unify the concepts of zero interaction in toxicology and statistics, 
Pharmacology and Toxicology , under review. 

5 

Loewe, S. and Muischnek, H. (1926) Uber kombinationswirkunger. I. Mitteilung: 
Hiltsmittel der gragstellung. Naunyn z Schmiedebergs, Arch. Pharmacol, 1 14, 313- 
326. 



10 Loewe, S. (1953) The problem of synergism and antagonism of combined drugs, 
Arzneimittle forshung. 3, 285-290. 

McCullagh, P. and Nelder, J.A. (1989) Generalized Linear Models (2nd ed.), New York: 
Chapman and Hall. 

15 

Meadows, S.L., Gennings, C, Carter, W.H., Jr., Bae, D.-S. (2002) Experimental Designs 
for Mixtures of Chemicals along Fixed Ratio Rays. Environmental Health Perspectives 
110(6):979-983. 

20 Moser, V.C. (1995) Comparisons of the acute effects of cholinesterase inhibitors using a 
neurobehavioral screening battery in rats. Neurotoxicol. Teratol . 17, 617-625. 

Murphy, S.D. and DuBois, K.P. Quantitative measurement of inhibition of the enzymatic 
detoxification of malathion by EPN (ethyl p-nitrophenyl thionobenzene phosphate). 
25 Proc. Soc. Exp. Biol. Med. 96:813-818, 1957. 

Nelder, J.A. and Mead, R. (1965). A Simplex Method for Function Minimization, The 
Computer Journal 7, 308-313. 

30 SAS Institute Inc, Version 8.2, Cary, NC. 

Teuschler, L., Klaunig, J., Carney, E., Chambers, J., Connolly, R., Gennings, C, Giesy, 
J., Hertzberg, R., Klaassen, C, Kodell, R., Paustenbach, D., Yang, R. (authors listed 
as Committee Chair, Co-Chair then alphabetically): Support of science-based 
35 decisions concerning the evaluation of the toxicology of mixtures: A new beginning. 
Regulatory Toxicology, 36:34-39, 2002. 

U.S. EPA (United States Environmental Protection Agency) (2000). Supplementary 
Guidance for Conducting Health Risk Assessment of Chemical Mixtures. EPA/630/R- 
40 00/002 



45 



89 



Appendix: 

Summary Statistics from Single Chemical and Mixture Experiments: 

The number of responders (r) represent the number of rats with gait scores >1 within each 
dose group. 



Group 


Dose 
(mg/kg) 


r 


N 


Single Chemical Data 


ACE 


0 


0 


8 


ACE 


3 


0 


8 


ACE 


10 


0 


8 


ACE 


30 


5 


8 


ACE 


60 


8 


8 


ACE 


120 


8 


8 










CPF 


0 


0 


20 


CPF 


1 


0 


12 


CPF 


2 


0 


8 1 


CPF 


5 


0 


12 


CPF 


10 


0 


20 


CPF 


20 


3 


8 


CPF 


25 


2 


12 


CPF 


30 


8 


8 


CPF 


50 


16 


20 










DIA 


0 


0 


16 


DIA 


5 


0 


16 


DIA 


25 


0 


16 


DIA 


50 


0 


8 


DIA 


75 


0 


16 


DIA 


125 


0 


8 


DIA 


150 


6 


8 


DIA 


250 


8 


8 










DIM 


0 


0 


8 


DIM 


5 


0 


8 


DIM 


10 


0 


8 


DIM 


25 


4 


8 


DIM 


50 


7 


8 


DIM 


75 


8 


8 










MAL 


0 


0 


8 


MAL 


100 


0 


8 


MAL 


500 


0 


8 
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Mixture data 



Group 


Total 
Dose 
(mg/kg) 


r 


N 


1 st Mixture Ray ('Full' Ray) a 


Positive Control data 


ACE 


6 


0 


6 


ACE 


30 


3 


14 


CPF 


10 


0 


6 


CPF 


20 


1 


8 


CPF 


30 


0 


6 


DIA 


50 


0 


6 


DIA 


125 


1 


14 


DIM 


10 


0 


6 


DIM 


25 


1 


8 


DIM 


50 


6 


6 


MAL 


350 


0 


14 


Mixture data 


CON 


0 


0 


14 


MIX 


10 


3 


12 


MIX 


55 


3 


12 


MIX 


100 


8 


12 


MIX 


200 


10 


12 


MIX 


300 


11 


12 


MIX 


450 


12 


12 


2 nd Mixture Ray ('Reduced' Ray) b 


Positive Control data 


ACE 


30 


3 


8 


CPF 


20 


0 


8 


DIA 


125 


0 


8 


DIM 


25 


3 


8 


Mixture data 


CON 


0 


0 


8 


MIX 


1.75 


0 


12 


MIX 


9.6 


2 


12 


MIX 


17.5 


2 


12 


MIX 


35 


6 


12 


MIX 


52.5 


12 


12 


MIX 


78.8 


11 


12 



a The total doses for the mixture studies are based on a fixed-ratio of the five pesticides. 
For the 'full' ray where all five pesticides are included the mixing ratio is (0.040: 0.002: 
5 0.03 1 : 0.825 : 0. 1 02) for (ACE: DIA: CPF: MAL: DIM), respectively. 

b For the 'reduced' ray the mixing ratio is (0.229: 0.01 1: 0.177: 0.000: 0.583) for (ACE: 
DIA: CPF: MAL: DIM), respectively. 
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Example 4. D s -OptimaI Designs for Testing Linear Hypotheses of Interest When Studying 
Combinations of Chemicals Using Multiple Fixed-Ratio Ray Experiments 
INTRODUCTION 

Organophosphorus (OP) pesticides are the most widely used pesticides in the Unites States 
5 (Aspelin, 1994). Colosio et al. (1999), Ma et al. (2002), and others have considered the effects 
of human exposure to pesticides. Due to the variety of agricultural, household, pet, and garden 
uses, there is high potential for multiple exposures from multiple routes. Thus, regulations that 
were based on the toxicity of these pesticides alone many not adequately protect humans from 
adverse health effects. The need for knowledge of how these compounds interact was 

1 0 highlighted in the 1 996 Food Quality Protection Act (FQPA), which directed the EPA to 
consider cumulative and aggregate exposures in the risk assessment process. 

The class of OP pesticides includes acephate, diazinon, chlorpyrifos, dimethoate, and 
malathion. These pesticides were chosen for study based on usage patterns (i.e., pesticides used 
on the same or similar crops) and market share (i.e., highest volume pesticides). Of the 

15 pesticides under study, malathion accounts for 82.5% of the relative dietary exposure to humans 
as projected by the U.S. EPA Dietary Exposure Evaluation Model (DEEM). In addition, 
preliminary analysis on available single chemical data indicated that malathion is not dose 
responsive (see Table 1) with respect to motor activity, a count of the number of passes made 
across a central area, which was the endpoint chosen to assess toxicity after exposure in rats. 

20 Since malathion does not appear to be dose responsive and it accounts for the largest portion of 
OP pesticides to which humans are exposed, it is of primary interest to study the effect of 
malathion on the remaining OP pesticides under consideration. Traditionally, to study 
interactions between malathion and the four active pesticides, researchers would consider 
factorial designs to study binary or tertiary chemical combinations. However, such studies that 

25 consider the effects of exposure to combinations of all five pesticides are not practical. 

As an alternative to studying chemical combinations using factorial designs, ray designs have 
been proposed to study polychemical mixtures. Brunden and Vidmar (1989) and others suggest 
the use of ray designs to support the estimation of a response surface. Alternatively, Gennings et 
al. (2002) and Meadows et al. (2002a) propose restricting inference to relevant fixed-ratio ray(s). 

30 Focusing inference to relevant fixed-ratio rays reduces the dimensionality and the amount of 
experimental effort associated with the study of chemical mixtures. Gennings et al. (2002) and 
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Meadows et al. (2002a) develop tests for departure from additivity using these fixed-ratio ray 
designs. Casey (2003) extended the work of Gennings et al. (2002) and Meadows et al. (2002a) 
to include hypotheses that test for interactions involving subsets of chemicals. Such hypotheses, 
described in Section 2, require the collection of data along a full fixed-ratio ray and a reduced ray 
5 where the chemical or subset of chemicals of interest (e.g. malathion) is removed from the 
mixture and the remaining chemicals are at the same relative ratios considered in the full 
mixture. 

The fixed-ratio rays used to study interactions between malathion and the four active 
pesticides were chosen based on relative dietary exposure estimates of each chemical as 

10 projected by the U.S. EPA Dietary Exposure Evaluation Model (DEEM). Consideration was 
given to a full ray given by (0.040: 0.002: 0.031: 0.102: 0.825) for (acephate: diazinon: 
chlorpyrifos: dimethoate: malathion) and a reduced ray given by (0.2286: 0.01 14: 0.1767: 
0.5833) for (acephate: diazinon: chlorpyrifos: dimethoate). Notice that the ratio of acephate to 
diazinon along the full fixed-ratio ray is (0.040: 0.002) which is a (20: 1) ratio. The ratio of 

15 acephate to diazinon along the second fixed-ratio ray is (0.2286: 0.01 14) which is also a (20: 1) 
ratio. Similar relationships hold for the other chemicals in the mixture; thus, the chemicals along 
the full and reduced rays are at the same relative ratios to one another. 

It is important to develop experimental designs for detecting and characterizing interactions 
among the pesticides. Assuming no interaction, it is of interest to determine dose locations and 

20 sample size allocation that provide precise parameter estimates and enough power to detect 
departure from additivity. In what follows we extend the work of Meadows (2001) and 
Meadows et al. (2002b) to develop designs for testing for the effect of subsets of chemicals on 
the mixture (i.e., to test for interactions involving subsets of chemicals) and for detecting 
departure from additivity across multiple fixed-ratio rays simultaneously. The criterion we 

25 consider aims to minimize the generalized variance of the parameters associated with the 

hypothesis of interest. While this criterion cannot claim to maximize the power of the hypothesis 
of interest, it should be related to an increase in power since it reduces the variance of the 
estimates involved. In order to determine experimental designs, appropriate models and 
hypotheses must be defined. The designs illustrated in this paper are based on the methods for 

30 detecting departure from additivity presented by Gennings et al. (2002), Meadows et al. (2002a), 
and Casey (2003) and are summarized in Section 2. 
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2. MODEL AND HYPOTHESES DEVELOPMENT 

Additivity Model: 

A generalized linear model is used to describe the relationship among the pesticides along the 
5 two fixed-ratio rays. In the quasi-likelihood framework, proposed by Wedderburn (1974), the 
generalized linear model relates the doses of the chemicals to the mean through a link function, 
g(m). Further it assumes that the variance of the response is of the form tV(m), where V(m) is a 
known function of m. 

Let K be the number of fixed-ratio rays under consideration for a mixture of c chemicals and 
10 a ( k) = [ai(k), «2(k), ,a C (k)], define the mixing ratio along the k th fixed-ratio ray, such that 

^ a i(k) = 1 . Let X(k)=[A:i(k), *2(k> *3(k), • • • -^c(k)] define a vector of doses at a given mixture point 

/=i 

along the k th ray and t =^ x i(k) define total dose. When ray designs are considered, total dose is 

the independent variable along the mixing rays and the amount of the i th compound in the 
mixture along the k th ray is given by a i(k )t. Following Meadows et al. (2002a), the additivity, i.e. 
15 no interaction, model along the k th fixed-ratio ray can be expressed as 
) = A> + M w ' + /V W + •- + Pc a c ( k)t 
= A>+(A«i<*) + 02*m + •••• + A<W (!) 

= a 

where: 

giMadd(kj) is the specified link function (see McCullagh and Nelder, 1989), 
/?o is the unknown parameter associated with the intercept, and 

20 9*(k) ~ X 0i a m ls me unknow 11 parameter associated with the slope along the k th ray. 

Interaction-Model-Based Method of Analysis: 

Meadows et al. (2002a) and Casey (2003) describe an interaction-model-based method of 
analysis where assumptions are made about the form of the underlying response surface. The 
25 significance of higher-order terms in the polynomial approximation of the dose-response 
relationship, expressed as a function of the total dose, indicates departure from additivity. 



94 



Meadows et al. (2002a) showed that with this approach, only mixture data along a fixed-ratio ray 
are necessary for detecting departure from additivity. Casey (2003) extended these results to 
detect departure from additivity across multiple fixed-ratio rays in the presence or absence of 
single chemical data. 

5 Using the quasi-likelihood framework, the interaction model along the k th ray can be 
expressed as, 

g(/W)) = a + o' m t + 0 2 y 2 + e; (k / + + (2) 

where: 

0'w =!>,«,■(*), 9l ( k) =£iiM(*) fl y(*)' t = Zli^V)V) fl w etc -> 

i=l 1=1 j=\ 1=1 j=l /=1 

i<j i<j<l 

10 g(jiimix(k)) is the specified link function (see McCullagh and Nelder, 1989), 

6 m is the unknown parameter associated with the first-order term, 

0* {k) for i = 2....c k is the unknown parameter associated with the i th way interaction, and 

Ck is the number of chemicals under study along the k th fixed-ratio ray. 

When single chemical data are available, the slope associated with each of the single chemicals 
1 5 are estimable. Thus, the model fit along the k th fixed-ratio ray is given by 

g(/W>) = A> + M (t )' + • • • + A* W + ^V 2 + 9 w/ + + • ( 3 ) 

The K fixed-ratio rays are generally fit simultaneously assuming a common intercept. If 
multiple control groups are experimentally evaluated, the assumption of a common intercept 
should be verified by comparing the means of the control groups. When single chemical data are 
20 available they are used in combination with mixture data to estimate the intercept and first-order 
terms and the higher-order terms are estimated using only the mixture data. In the absence of 
single chemical data, all of the model parameters are estimated using the mixture data. 
Parameter estimates are found by maximizing the quasi-likelihood (e.g., McCullagh and Nelder, 
1989). 

25 Under the hypothesis of additivity, the parameters associated with interaction along the K 

fixed-ratio rays, namely (0 2 * (I) , <9 3 * (I) , 9* 2(2) , 0 3 * (2) <9 C * 2(2) > ,0l w > e l(K) >~fl K m ) > 816 

zero. Define the pxl vector = [/? 0 , 0* (1)J 0 2 * (1) ,.•■> 0* Ci(l) ,....,9l {K) ,0* 2(K) , 6* Ck(K) ]' as a vector of 
model parameters and define b a dd as a matrix of contrasts (zeros and ones) such that 
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An overall hypothesis of additivity, based on the significance of higher-order terms, is given by 
H„:IWY=0 (4) 
In addition to the hypothesis described above, Casey (2003) developed hypotheses for 
5 detecting and characterizing interactions due to subsets of chemicals. The hypothesis that a 

chemical or subset of chemicals is not involved in interactions with the remaining components of 
the mixture is given by 



2(/«/7) °2(reduced) 



(«!(>«) ) 2 


(° '{(reduced)) 




^i(reduced) 


Kmo) 3 


( a \(reduced)) 




@m(reduced) 



( a i(/u//)) ( a i(reduced)) 



where the reduced ray (parameters denoted with subscript 'reduced') represents the mixture 
1 0 study where the chemical or subset of chemicals of interest is removed from the mixture and the 
remaining components are studied at the same relative ratios as given in the full ray. Parameters 
along the full ray, where all of the chemicals are considered in the mixture, are denoted with 
subscript 'full'. This hypothesis compares interaction parameters along the reduced ray (e.g. the 
no malathion ray) to interaction parameters along the full ray. 
1 5 For example, consider a mixture study with three chemicals. Let (fly^u//;: a 2 (f U uf. ci3(fuii)) denote 
the mixing ratio. The model along the fixed-ratio ray is given by, 

giMifaii)) = Po +0 Kjuin t + 02(/«//)' 2 +^3(>//) f3 

where: 

^l(fall) = P\ a \(full) + Pi a 2(full) + /?3 a 3(>//)» 

02(/ull) = Pl2 a \(full) a 2{full) + Pu a i(fu!l) a l(futl) + P 2i a 2( Jull^if"'') ' ^ 
^i(fiill) ~ P\23 a \(Jutl) a 2(faU) a 3(fall)- 

20 Suppose the third-order term is not significant and it is of interest to determine the effect of the 
third chemical on the mixture. Thus, data are collected along a second fixed-ratio ray given by 
(a t(re d U ced)'- a 2 ( r educed)) where the third chemical is removed from the mixture and 
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= induced) ^ m()del fit along thig reduced fay ig giyen 

a Kfull) a j(reduced) 



git* (reduced)) ~ flo + ^(reduced)* + ^(reduced/ 



where: 



^l(reduced) = 
@l{reduced) ~ Pu a \(reduced) a 2(reduced)- 

Since the reduced experiment is performed at the same relative ratios as those in the original or 
mil ray, 

el 



u l(reduced) u Unreduced) _ u 2(reduced) 
( a \(full)) ( a \(reduced)) ( fl 1( reduced ) ) 

2 becomes 

K/wo) 



6 

and — 2(/ "'° , becomes 



/ x2 ~~ / 2 x2 , x2 ^^23", . 2 

\ a \(full)) \ a \{reduced)) \. a \(full)) K a \( full)) 

f2(jm Onreduced) _ R Q XJulD^iM) a a 2(full) a 3(full) 

( \ 2 ( \ 2 ~ 13 ( \ 2 ^ 23 ( \2 ' 

\ a i(fiill)) \ a i(reduced)) \ a i(full)) \ a i(full)) 



The hypothesis, given by H 0 



= 0, tests for the effect of the third 



_( a i(/uH)) ( a i{reduced)) 

chemical on the mixture (i.e. the significance of /? u and/or 0 2 3)- If there is sufficient evidence to 
reject this hypothesis, we conclude that the third chemical is involved in at least one two-way 
interaction. Further experimentation would be necessary to determine specifically which 
chemicals are involved in these interactions. If we fail to reject the hypothesis and the study is 
reasonably powered, we can conclude that the third chemical does not interact with the 
remaining components of the mixture. 

A Wald-type test for testing the hypotheses given in (4) and (5), is given by 

w= (b7) , [b«bT'(b7) ) (6) 
Mt 
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where zB is the variance-covariance matrix of 7 and b is any contrast matrix (e.g. b a dd; bj„teract)- 
Since 7 is distributed asymptotically normal (McCullagh and Nelder, 1989) with mean y and 
variance zfl , it follows that W is approximately distributed chi-square with M degrees of 
freedom (M is the number of tests considered simultaneously or the number of rows in b). The 
moment estimate for t is expressed as 



1 



x 2 



where X 2 is the generalized Pearson statistic (McCullagh and Nelder, 1989), which is 
asymptotically distributed chi-square with N-p degrees of freedom. In the quasi-likelihood 
framework, McCullagh (1983) defines the large sample variance-covariance matrix for 7 as 
r[/(7)] _1 , where 1(g) is the expected quasi-information matrix. Let N be the number of 
observations. McCullagh (1983) expressed the expected quasi-information as 
/(7) = D'OWD 

where 



(7) 



(8) 



and V(//) = 



VM 0 
0 VQi 3 ) 



0 



0 



0 
0 
0 



15 Replacing Twith y in (8), Q = [I(y)] 1 is a consistent estimate for W (McCullagh, 1983). 

Replacing t with f and W with Q in (6), W is approximately distributed F with M numerator 
degrees of freedom and N-p denominator degrees of freedom, where p is the number of 
parameters. 

3. DATA DESCRIPTION 

20 The five-pesticide study described in the introduction is used to illustrate the design 
methodology developed in the following sections. Due to the high potential for aggregate 
exposure to OP pesticides, researchers were initially interested in testing the hypothesis of 
additivity, given in (4), along the full five-pesticide fixed-ratio ray, given by (0.040: 0.002: 
0.031: 0.102: 0.825) for (acephate: diazinon: chlorpyrifos: dimethoate: malathion). Single 

25 chemical data were collected and used to estimate an additivity model along the ray (results 

shown in Table 1). Since motor activity response is a count variable, it was assumed that V(/x)= 
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(i and the log link function (McCullagh and Nelder, 1989) was used to fit the generalized linear 
additivity model given in (1). A plot of this additivity model, seen in Figure 10, was useful in 
specifying the active dose range of 0 mg/kg to 450 mg/kg along the ray. Prior to the design 
methods developed here, an equally spaced design with equal allocation and extra dose points in 
the low dose region was chosen to study interactions along the full ray. Eighty-six animals were 
considered. Fourteen animals were evaluated in the control group. Twelve animals each were 
evaluated at 10, 55, 100, 200, 300, and 450 mg/kg. Summary statistics for the single chemical 
and mixture data are provided in the appendix. 

Since the mixture contains five chemicals the single chemical and mixture data were fit to the 
fifth-order generalized linear model given in (3). Maximum quasi-likelihood estimates were 
found using PROC GENMOD with DIST=POISSON and LINK=LOG in SAS® (version 8.2) 
and are provided in Table 1 . A plot of the interaction model is provided in Figure 11. The 
hypothesis given by 



20 



H 0 :b add 7= 



y 2<Jull) 
7 3(futt) 



= 0, 



(9) 



where 7 = 


[A 










"0 


0 


0 


0 


0 0 10 0 




0 


0 


0 


0 


0 0 0 1 0 


b add = 


0 


0 


0 


0 


0 0 0 0 1 




0 


0 


0 


0 


0 0 0 0 0 



, was used to test the hypothesis of additivity. With a 



p-value of 0.007, there was sufficient evidence to reject this hypothesis and conclude departure 
from additivity exists along the full five pesticide fixed-ratio ray. The fifth-order interaction 
term was not significant (p-value, 0.22); thus, it was removed from the model. In addition, the 
parameter associated with malathion was removed from the model due to lack of significance (p- 
value, 0.85) which indicates the lack of activity of malathion. The fourth-order interaction term 
was marginally significant with a p-value of 0.07. Although marginal, we may conclude that up 
to fourth-order interactions exist along the full ray; however, without further experimentation, we 
cannot determine which chemicals are involved in these interactions. 
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Table 1. Estimated model parameters where the single chemical and mixture data along the mil 
fixed-ratio ray are fit to the model given in (3) with the log link (allowing V(m)=m). 



Parameter 


Estimate 


SE 


p-value 


Po 


5.326 


0.020 


<0.0001 


Pi (Acephate) 


-0.018 


0.001 


O.0001 


p 2 (Diazinon) 


-0.004 


0.0004 


O.0001 


P 3 (Chlorpyrifos) 


-0.025 


0.002 


<0.0001 


p 4 (Dimethoate) 


-0.0123 


0.001 


<0.0001 


0* (2 nd Order Interactions) 


-3E-5 


1E-5 


0.0097 


0\ (3 rd Order Interactions) 


1.5E-7 


7.2E-8 


0.0345 


B\ (4 th Order Interactions) 


-1.9E-10 


1E-10 


0.0714 


X 


3.503 







NOTE: ^ fii a i(f U ii) = 0'(faio = -0.0028. The fifth order interaction term was removed from the 

5 model due to lack of significance (p-value, 0.2203). Similarly, the slope for malathion was 
removed from the model (p-value, 0.8485). 

The existence of fourth-order interactions along the full fixed-ratio ray and the lack of activity 
of malathion provided motivation for studying the effect of malathion on the four active 

10 pesticides. That is, it is of interest to determine if malathion interacts with the other four 

pesticides even though it is not active alone. This can be achieved through studying a reduced 
ray, where the mixing ratios are given by (0.2286: 0.0114: 0.1767: 0.5833) for acephate, 
diazinon, chlorpyrifos, and dimethoate. In this mixture, malathion was and the remaining 
pesticides are at the same relative ratios as in the full ray. Thus, it is of interest to determine a 

1 5 design along the reduced ray associated with the hypothesis that malathion does not interact with 
the four active pesticides. This hypothesis is given by 



r 2(full) °2(reduced) 



Kjwo ) 2 


( a \(reduced)) 


®Hfull) 


Unreduced) 




( a i(reduced)) 




@4(reduced) 



( a i(full)) ( a i(reduced)) 



where y = [fi 0 , /?, , fi 2 , , p A , 6* 2(m , 0\ (m , 6\ (full) , 0 2 \ reduced) , 0* 3(reduced) , 0l (reduced) ] 'and 
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0 0 0 0 0 — 0 0 i 0 0 

(0.04) 2 (0.2286) 2 

0 0 0 0 0 0 — l — 0 0 0 

(0.04) 3 (0.2286) 3 

0 0 0 0 0 0 0 — l — 0 0 



(0.04) 4 (0.2286) 4 

The objective of this paper is to provide methodology to determine at what dose levels and in 
what proportion of the total sample size observations should be taken. 
4. EXPERIMENTAL DESIGNS 
5 D s -optimality is a special case of D-optimality, defined by Kiefer and Wolfowitz (1959). D- 
optimal designs minimize the generalized variance of the model parameters, where the 
generalized variance is defined as the determinant of the variance-covariance matrix. While this 
criterion is not directly related to optimal hypothesis test power, use of such designs should result 
in increased power since the parameters of interest are estimated more precisely. D s -optimality 
10 minimizes the generalized variance for a subset of model parameters. We have seen in (4) and 
(5) that certain hypotheses of interest in chemical mixture problems can be expressed in terms of 
subsets of the parameter vector. D s -optimality is appropriate when the precision of a subset of 
parameters is of interest (Atkinson and Donev, 1992). 

In the quasi-likelihood framework, the determinant of the variance-covariance matrix is given 
15 by jzft -1 1 , where W is the expected quasi-information matrix given in (8). When the hypothesis 

of interest is given by a subset of parameters (e.g. the hypothesis of additivity given in (4)) or is 
based on linear combinations of model parameters (e.g. the hypothesis given in (5)), the variance 
covariance matrix is given by r(bfl b') -1 , where b is an appropriate matrix of contrasts (e.g. b ad d; 
binteract). The expected quasi-information matrix depends on parameter values, total dose values, 
20 and the allocation of observation to the dose groups. Given expected parameter values, the total 
sample size, and the form of the hypothesis of interest, the Nelder Mead algorithm (Nelder and 
Mead, 1965) can be used to determine the total dose levels and the allocation of observations that 
minimizes the generalized variance, |r(bfl b')" 1 1 , for the model parameters of interest. 

Using the methods described above, it is of interest to determine the design along the reduced 
25 ray which minimizes the generalized variance associated with the parameters involved in the 
hypothesis, given in (10), that malathion dose not interact with the four active pesticides. The 
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D s -optimal design methods condition on the total sample size, the model parameters, and the 
number of dose groups of interest. Thus, these values must be specified prior to implementing 
the methods developed here. 

Values for the intercept (/So), single chemical slopes (/? i,02,03, fid, interaction parameters 
5 along the full ray ( 9[ {fuU) , 9* XfilU) , 0\ {full) ), and the variation parameter (t) are provided from the 
analysis performed on the full ray (Table 1). Preliminary mixture data along the reduced fixed- 
ratio ray would be ideal for specifying the additional higher-order interaction terms 
( Reduced) > Reduced) > Reduced) )• However, since these data are not available the higher-order 
interaction parameters along the full ray are used as a guide to specify similar parameters along 
1 0 the reduced ray under the alterative hypothesis that malathion is involved in interactions with the 
four active pesticides. Since 0* 2(filU) , 0* Xfull) , Q\ m ,ai<fuii), and aj( re d uce d), have been specified we 
can define a model along the reduced ray under the null hypothesis (i.e., under the assumption 
that malathion does not interact with the remaining pesticides). For example, consider the case 
where malathion is not involved in any three-way interactions. Under this assumption 



( a \{full)) (- a \(reduced)) 



(11) 



( a \(full)) ( a i{reduced)) 

^ 3 </""> *(g 
/ \3 \ u l{reduced)J ) v ^(reduced) 

\ a \(full) ) 

Parameters values along the reduced ray under the assumption that malathion does not interact 
with the active pesticides are provided in Table 2. 

Table 2. Model parameters along the reduced fixed-ratio ray under the null hypothesis (i.e. 
20 malathion does not interact with the remaining pesticides) and the alternative hypotheses that (a) 
malathion is involved in two-way interactions, (b) malathion is involved in two and three-way 
interactions, and (c) malathion is involved in two, three, and four-way interactions. 

Parameter No Malathion Two-way Two and Two, Three and 
Interactions Only Three-way Four-way 

|3o 5.326 5.326 5.326 5.326 

0i(Acephate) -0018 -0.018 -0.018 -0.018 

0 2 (Diazinon) -0.004 -0.004 -0.004 -0.004 

0 3 (Chlorpyrifos) -0.025 -0.025 -0.025 -0.025 

04 (Dimethoate) -0.012 -0.012 -0.012 -0.012 
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0\ (2 nd Order Interactions) 


-9.8E-4 


-9.5E-4 


-7.84E-4 


-6.664E-4 


#3 (3 rd Order Interactions) 


2.8E-5 


2.8E-5 


2.52E-5 


1.96E-5 


9\ (4 th Order Interactions) 


-2.0E-7 


-2.0E-7 


-2.0E-7 


-1.4E-7 


t 3.503 


NOTE: £M (reM = 


°- 01 









Using the no malathion interaction case as a guide, we can specify values of model parameters 
under the alternative hypothesis based on cases that are of interest to investigators. First, we 
5 consider the case that malathion is involved in two-way interactions. In this case the third and 
fourth-order interaction terms are determined following the process defined in (1 1) and the 
second-order term is defined based on changes in the mean values that researchers define as 
biologically meaningful. We also consider the case where malathion is involved in two and 
three-way interactions and the case where malathion is involved in two, three, and four-way 
10 interactions. The interaction model parameters along the reduced ray for these three cases are 
provided in Table 2. 

In addition to using the parameter values under the null hypothesis or no malathion interaction 
case, plots of the interaction models under the alternative hypothesis and tables of total dose 
values where a given percent change in the mean value occurs between the no malathion 

1 5 interaction case and the specified alternative cases are useful in defining biologically meaningful 
interactions. The plots are provided in Figures 2(a-c). The additivity model along the reduced 
ray is also provided to aid in the interpretation of the interaction model specified under the 
alternative hypothesis. Table 3 provides the doses where a given percent change in the mean 
response is observed between the no malathion interaction case and the model specified under 

20 the alternative hypothesis. Since the generalized linear model is nonlinear, the differences in the 
curves are not constant. For example, if we are interested in detecting at least a 5% change in the 
mean for the two-way interaction case, the minimum total dose value that is associated with such 
a change in the mean is 40.33 mg/kg. 
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Table 3. Lowest total dose values along the reduced fixed-ratio ray that result in 5%, 
10%, 15%, and 20% changes in mean responses between the model under the null 
hypothesis and the model under the alternative hypotheses that (a) malathion is involved 
5 in 2-way interactions, (b) malathion is involved in 2 and 3-way interactions, and (c) 
malathion is involved in 2, 3, and 4-way interactions. 



Change in 


Two-way 


Two and Three-Way 


Two, Three, and 


Mean 


Interactions 


Interactions 


Four- Way 








Interactions 


5% 


40.33 


18.37 


15.78 


10% 


56.36 


28.71 


27.04 


15% 


68.25 


42.98 


77.69 


20% 


77.96 


81.88 


81.46 



Thus far, we have used available information provided by the analysis on single 

10 chemical and mixture data along the full ray, in combination with plots, to specify model 
parameters associated with the alternative hypotheses of interest. Now we must specify 
the size of the design and the available sample size. Suppose researchers indicate that a 
maximum sample size of 80 animals is possible for the mixture study along the reduced 
ray. For a mixture of c chemicals, (c+1) is the minimum number of dose group necessary 

1 5 to detect departure from additivity along a fixed-ratio ray. Thus, the minimum design 
along the reduced ray is a five-point design. 

The Nelder-Mead algorithm (Nelder and Mead, 1965) is used to determine the optimal 
allocation of the 80 available animals to five total dose groups that minimizes the 
generalized variance associated with the parameter values related to the hypothesis that 

20 malathion is not involved in interactions with the remaining four pesticides. The Nelder- 
Mead algorithm requires starting values for the dose levels and the sample size 
allocations. For each of the three cases considered we force the design to include a 
control group and a range of starting values (0 to 80 by 5) for the remaining total dose 
values were considered. Initially, equal allocation across the dose groups was assumed. 

25 The D s -optimal designs are provided in Table 4. Notice that the design for the case 

where it is assumed that malathion is involved in two-way interactions is the same as the 
design for the case where malathion is assumed to be involved in two, three, and four- 
way interactions. Also, notice that in each of the three cases the allocation of the total 
sample is relatively equal across the dose groups. Finally, if we look at the plots of the 
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models under the alternative hypotheses, provided in Figures 2(a-c), we see that the total 
dose levels determined by the D s -optimal designs are located where the curvature exists 
along the model considered under the alternative hypothesis and where the largest 
differences occur between the no malathion interaction model and the model under the 
5 alternative hypothesis. 

Table 4. Five point D s -optimal designs for the hypothesis, given in (10), that malathion 
does not interact with the four active pesticides. The top row corresponds to the total 
dose locations and the bottom row of the design corresponds to the sample size 
10 allocation. Parameter estimates along the full fixed-ratio ray are provided in Table 1. 
Assumed model parameters under the alternative hypothesis along the reduced ray for 
each of the three cases considered are provided in Table 2. 





D s -Optimal Design 


Two-way Interaction Case 


-I 


0 33 35 65 801 
16 16 16 16 16J 


Two and Three-way Interaction Case 


H 


r 0 30 35 61 801 
16 16 16 16 16j 


Two, Three, and Four-way Interaction 
Case 


H 


0 33 35 65 801 
16 16 16 16 16J 



1 5 For each of the three cases considered above, we found the D s -optimal design with the 
minimum number of design points, c+1 . It is of interest to examine the effect of adding 
additional dose groups. For example, consider the six, seven and ten point D s -optimal 
designs, provided in Table 5, for the alternative hypothesis that malathion is involved in 
two-way interactions only. Recall that each of the designs forces a control group. Notice 

20 that the five-point design places two dose groups in the middle of the dose region at the 
first curve in the model and two doses in the high dose region where there is additional 
curvature and the largest difference between the null and alternative models occurs. The 
six-point design continues to follow this pattern. In fact the six-point design is not much 
different from the five-point design as 79 mg/kg is not much different from 80 mg/kg. 

25 As more design points are added, the doses are added to the middle and high dose regions 
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and are relatively equally spaced in these dose regions. Figure 12 illustrates the 
placement of dose points for the five, six, seven, and ten-point designs given in Table 5. 



Table 5. Six, seven, and ten point D s -optimal designs for the hypothesis, given in (10), 
that malathion does not interact with the four active pesticides and the alternative 
hypothesis that malathion is involved in two-way interactions. The top row of the design 
corresponds to the total dose locations and the bottom row corresponds to the sample size 
allocation. Parameter estimates along the full fixed-ratio ray are provided in Table 1 . 
Assumed model parameters under the alternative hypothesis along the reduced ray are 
provided in Table 2. 



Two-Way Interaction 

Case 


D s -Optimal Design 


Six Point Design 


JO 30 35 60 79 80l 
[13 13 13 14 13 13J 


Seven Point Design 


JO 25 31 36 56 66 80| 
[10 12 12 12 12 12 12J 


Ten Point Design 


JO 26 30 35 40 55 60 65 79 80] 
[8 888888888J 



5. EFFICIENCY 

When it is of interest to detect and characterize interactions among chemicals in a 
mixture, multiple hypotheses are often considered. Different hypotheses can result in 
different designs. In choosing the appropriate design, consideration can be given to the 
primary hypothesis of interest. Given the design corresponding to the primary hypothesis 
of interest (*F) and the design corresponding to the secondary hypothesis OF*), efficiency, 
described by Atkinson and Donev (1992) and given by 



Eff = 



\M{W)\ 



\M{ ¥ ') 



can be used to determine how efficient T is relative to *F for testing the secondary 
hypotheses (p is the number of model parameters). For example, we are primarily 
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interested in testing for interactions between malathion and the four active pesticides. 
However, in order to ensure appropriate conclusions are made with respect to the effect 
of malathion, it is important to consider the hypothesis for additivity, given in (4), along 
the reduced ray. If the design chosen to study interactions between malathion and the 

5 remaining components of the mixture does not adequately allow for the detection of 
higher-order terms along the reduced ray, conclusions about the effect of malathion may 
be misleading. Thus, it is important to consider the efficiency of the designs given in 
Table 4 for testing the hypothesis of additivity given in (4). 

Let \|/interact represent the design associated with the hypothesis, given in (10), that 

10 malathion is not involved in any interactions and let i|/ a dd represent the design for the 
overall hypothesis of additivity along the reduced ray, given in (4). Following Atkinson 
and Donev (1992), the D s -efficiency of i|/i n teract for testing the hypothesis of additivity 
along the reduced ray can be expressed as 



the hypothesis of additivity when the v|/i n teract design is considered, |M(^ add )| is the 
determinant of the variance-covariance matrix associated with the hypothesis of 
additivity when the vy a dd design is considered, and p is the number of parameters. If D s . e ff 
equals one, ^interact is fully efficient for testing hypothesis of additivity. The smaller D s . e ff 

20 the less efficient interact is for testing the hypothesis of additivity. 

Consider the D s -optimal design (^interact), given in Table 4, for the case where it is 
assumed that malathion is involved in two, three, and four-way interactions. We are 
interested in determining the efficiency of this design for testing the overall hypothesis of 
additivity. The D s -optimal design for testing for departure from additivity along the 

25 reduced ray is given by 



The efficiency of v|/ interact for testing for the presence of the higher-order terms along the 
reduced ray is given by 




(12) 



15 where \M (^ interact )J is the determinant of the variance-covariance matrix associated with 




(16 16 16 16 16 



0 30 35 65 80^ 
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Thus, the D s -optimal design associated with testing interactions between malathion and 
the four active pesticides under the assumption that malathion is involved in two, three, 
and four-way interactions is efficient for testing for the presence of higher-order terms 
5 along the reduced ray. Notice that these two designs are similar to one another, thus we 
would expect high efficiency. 
6. DISCUSSION 

The goal of this paper was to develop methodology for determining dose locations and 
sample size allocation that provide precise parameter estimation and enough power to 
10 detect departure from additivity and interactions due to subsets of chemicals along fixed- 
ratio rays. The precision of parameter estimates is an important issue as it is associated 
with increased power. Consideration has been given to D s -optimal designs which are 
concerned with minimizing the generalized variance associated with the hypothesis of 
interest. 

1 5 The hypotheses developed in Section 2 require the fit of higher-order polynomial 
models. Working with such models can cause numerical problems, particularly as the 
order of the polynomial increases. For example consider the model parameters that 
results from fitting the generalized linear model with the log link to the single chemical 
data and to the data along the full fixed-ratio ray (estimates provided in Table 1). Notice 

20 that the parameter associated with the fourth-order term is -1 .9E- 1 0 with a standard error 
of 1.0E-10. We are considering a dose range of 0 mg/kg to 450 mg/kg along the full ray 
and 0 mg/kg to 80 mg/kg along the reduced ray. Raising such dose values to powers of 
greater than or equal to two results in relatively large numbers; thus, parameters estimates 
associated with higher-order terms become smaller as the order of the polynomial term 

25 increases. In some situations such small numbers cause calculation problems. It should 
be noted that the degree of the polynomial model chosen follows directly from the 
number of components in the mixture (Meadows et al., 2002a; Casey (2003)). Therefore, 
these higher-order polynomial terms are important in detecting and characterizing 
interactions. One proposed solution to ease numerical problems is to consider scaling the 

30 dose range. For example we could divide the doses by 10 or 100 and work with either 
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centigrams or decigrams instead of milligrams. Although the results are not presented 
here we examined the effect of scaling the dose on the D s -optimal designs. For the 
examples we considered the designs reported here were equivalent to the designs along 
the scaled dose region. 

5 Since the hypotheses of interest involve only a subset of the model parameters, Ds- 
optimal designs were developed for fixed-ratio rays. The D s -optimal examples presented 
here considered the design along a single fixed-ray ray. These methods can be extended 
to simultaneously consider dose location and sample size allocation for multiple rays. 
Although the examples presented here resulted in relatively equal allocation it is 

1 0 important to note that equal allocation is not always optimal. D s -optimal designs 

condition on model parameter values specified under the alternative hypothesis. In the 
case where competing hypotheses are of interest, methods for determine the D s -efficiency 
of one design relative to another were also developed. 

The designs presented here were developed for the generalized linear model. In 

15 evaluating the risk associated with exposure to mixtures it is often of interest to detect a 
threshold (Schwartz et al., 1995). The methods for detecting overall departure from 
additivity and for detecting interactions between subsets of chemicals and the remaining 
components of a mixture are readily extended to the threshold model (Casey, 2003). The 
methodology for determining D s -optimal designs can readily be applied when threshold 

20 models are of interest. 
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APPENDIX 



Summary statistics for motor activity response associated with the single chemical data. 



Chemical 


Dose 


Mean Motor 


Standard 


Sample Size 




mg/kg 


Activity 


Deviation 








Response 






Acephate 


0 


217.88 


35.77 


8 




3 


200.13 


34.13 


8 




10 


165.88 


25.02 


8 




30 


108.75 


62.51 


8 




60 


58.25 


24.23 


8 




120 


33.25 


27.41 


8 


Diazinon 


0 


206.69 


34.77 


16 




5 


190.88 


28.49 


16 




25 


215.56 


24.02 


16 




50 


183.13 


24.44 


8 




75 


165.69 


33.05 


16 




125 


152.00 


38.65 


8 




150 


76.5 


35.4 


8 




250 


61.25 


47.07 


8 


Chlorpyrifos 


0 


190.00 


14.36 


8 




2 


192.75 


38.61 


8 




10 


172.13 


17.55 


8 




20 


157.13 


28.31 


8 




30 


80.13 


34.85 


8 




50 


56.38 


29.98 


8 


Malathion 


0 


195.86 


19.06 


7 




100 


201.5 


28.38 


8 




500 


203.75 


28.34 


8 


Dimethoate 


0 


195.75 


33.12 


8 




5 


188.25 


24.21 


8 




10 


188.63 


53.05 


8 




25 


107.75 


37.23 


8 




50 


103.75 


51.59 


8 




75 


101.50 


59.57 


8 
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5 

Summary statistics for motor activity response associated with mixture data along the full 



fixed-ratio ray. 

Total Mean Motor Standard Sample Size 
Concentration Activity Deviation 
Dose (mg/kg) Response 

0 " 199.43 " " 20.77 14 
10 200.92 27.10 12 
55 167.92 37.55 12 
100 117.08 47.34 12 
200 95.17 34.03 12 
300 72.25 40.93 12 
450 6O08 46J6 12 
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Example 5. Statistical Analysis of Interactive Cytotoxicity in Human Epidermal Keratinocytes 

Following Exposure to a Mixture of Four Metals 

INTRODUCTION 

Both occupational and environmental exposures to hazardous metals are significant 
toxicological concerns. Not only do metals such as arsenic (As), cadmium (Cd), chromium (Cr), 
lead (Pb) and nickel (Ni) lead to toxicity in exposed individuals, but epidemiological evidence 
suggests that some, if not all, of these metals are human carcinogens. Many laboratories, using a 
variety of experimental systems and techniques, have invested substantial resources into 
mechanistic studies of metal toxicity. We are, however, still a long way from understanding the 
relationship between these observed cellular effects and the toxicity of metals in exposed 
populations. To further complicate the picture, individuals are rarely exposed to single metals. 
Rather, they are exposed to metals in combination with one another or with other environmental 
contaminants such as those prevalent at hazardous waste sites. In fact, 100 or more different 
chemicals can be found at a single hazardous waste site in varying combinations in water, soil, 
and air (De Rosa et al, 1996). Like most chemical mixtures, however, the toxicity of these metal- 
containing mixtures has not yet been studied. It is anticipated that many of the toxic effects of 
metals, including carcinogenicity, can be modified by concurrent exposure to other metals. 
We are interested in understanding the biological processes that drive metal-mediated 
cytotoxicity and carcinogenesis, particularly with respect to metals in environmentally relevant 
chemical mixtures. The metals chosen for these studies were As, Cd, Cr, and Pb, which are 
among the top contaminants in site frequency count by the ATSDR Completed Exposure 
Pathway Site Count Report (ATSDR, 1997); three of these, As, Pb, and Cd are among the 
Superfund Top 10 Priority Hazardous Substances (De Rosa et al, 1996), i.e. those considered to 
pose the greatest hazard to human health. In addition, as demonstrated by ATSDR using the 
HazDat database, these metals are present in eight (five) of the Top 10 Binary Combinations of 
Contaminants in soil (water) (Fay and Mumtaz, 1996). Mechanistic studies have demonstrated 
that these four metals affect multiple intracellular targets, (including both nucleic acids and 
proteins) and exert a variety of diverse effects on cells in culture (including DNA damage and 
mutagenesis, enzyme inactivation and altered metabolism, and aberrant signal transduction). The 
role played by these intracellular alterations in either metal-mediated acute toxicity or 
carcinogenicity is currently uncertain. 
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Epidermal keratinocytes, both in culture and in vivo, have been used extensively to study 
the effects of chemical exposure at the cellular level, particularly with regards to malignant 
transformation. These types of studies have defined phenotypic and genotypic alterations that 
occur during transformation of normal skin cells (Pietenpol, Holt, Stein, and Moses, 1990; 
Gaido, Maness, Leonard, and Greenlee, 1992; Glick, Sporn, and Yuspa, 1991; Glick, Lee, 
Darwiche, Kulkarni, et al, 1994; Choi, Toscano, Ryan, Riedel, et al, 1991; Punnonen, Denning, 
Rhee, and Yuspa, 1994). In addition, for As at least, the skin is a critical target organ for metal- 
mediated carcinogenesis and other proliferative disorders such as hyperkeratosis. Due to high As 
concentrations in many drinking water supplies, this has become a health problem of global 
proportions. Arsenic has substantial effects on epidermal keratinocytes in vitro, altering 
expression of several growth regulatory factors such as TGFa and TGFb, and inhibiting the 
normal process of differentiation (Germolec, Yoshida, Gaido, Wilmer, et al, 1996; Yen, Chiang, 
Chang, Tsai, et al, 1996; Germolec, Spalding, Boorman, Wilmer. et al, 1997; Kachinskas, Qin, 
Phillips, and Rice, 1997). This latter effect is also seen in keratinocytes with Cr administration, a 
well known skin sensitizer (Cohen, Kargacin, Klein, and Costa, 1993; Ye, Zhang, Young, Mao, 
and Shi, 1995). We have, therefore, chosen human keratinocytes for our toxicological studies of 
metal mixtures. These studies have two distinct phases and are being carried out in both primary 
human epidermal keratinocytes and virally or spontaneously immortalized keratinocyte cell 
lines: 1) evaluation of the acute cytotoxicity of As, Cd, Cr, and Pb in this cell type, both 
individually and in combination in a four-metal mixture, and 2) determination of the 
carcinogenic potential of the metals and metal mixture in low-dose long-term chronic exposure 
scenarios. In both phases of our work a great deal of emphasis is being placed on identifying 
toxicological interactions among the metals, whether synergistic or antagonistic. Subsequent 
mechanistic studies will increase our understanding of these interactions and aid in development 
of more realistic risk assessment strategies for metal-containing chemical mixtures. As part of 
this larger study, we present here a detailed statistical analysis of the cytotoxic interactions of As, 
Cd, Cr, and Pb in normal primary human keratinocytes. 

Of primary importance in studying these four metals is the determination and 
characterization of their interactions when present in a mixture. For example, it is of interest to 
determine if the metals interact in such a way as to increase the toxic effect above what one 
would expect to observe from any single chemical. If the metals do not interact, it is said that 
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they may either be functionally independent or may act in a dose-additive fashion. One 
definition of additivity is given by Berenbaum (e.g., 1985) and is based on the classical 
isobolograms for the combination of two chemicals (e.g., Loewe and Muischnek, 1926; Loewe, 
1953). That is, in a combination of c chemicals, let Xj represent the concentration/dose of the i tf 
component alone that yields a fixed response, y 0 , and let xj represent the concentration/dose of 
the i th component in combination with the c agents that yields the same response. According to 
this definition of additivity if the substances combine with zero interaction, then 




The additivity surface can be used to test the null hypothesis that the chemicals act in an 
additive fashion. The experimental data necessary and sufficient to support the estimation of the 
additivity surface are single chemical dose-response data (e.g., Gennings, et al., 1997). The 
experimental design considered herein also includes a fixed ratio of the four metals at varying 
total dose levels. It is of interest to determine if the observed mean responses from these mixture 
points coincide with the predicted response under the assumption of additivity. This general 
approach is described in detail by many authors, including Berenbaum (1985), Kelly and Rice 
(1990), Gennings and Carter (1995), Gennings et al (1997). However, these authors do not 
address the issue of simultaneous testing at multiple mixture points. 

In what follows, two simultaneous tests for departure from additivity are developed. If 
the overall test of additivity is rejected, then corrected single degree of freedom tests can be used 
to determine at which mixture point(s) an interaction exists. The correction is based on 
Hochberg's (1990) sequentially rejective approach. Prediction intervals are then used to 
characterize the interaction as synergistic or antagonistic. In addition, a simultaneous confidence 
band on the difference in the predicted response and that predicted under additivity is compared 
to a zero reference level to identify concentration ranges of departure from additivity. 
MATERIALS AND METHODS 

Chemicals. Sodium metaarsenite (NaAs02), cadmium chloride (CdCh), chromium oxide 
(Cr0 3 ), chromium chloride (CrCl 3 ), lead acetate [(C 2 H 3 02)2Pb 3H 2 0], and 3-[4,5- 
Dimethylthiazol-2-yl]-2,5-diphenyltetrazolium bromide (MTT) were purchased from Sigma 
Chemical Co. (St. Louis, MO). 
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Cells and culture conditions. Cryopreserved, pooled normal human epidermal keratinocytes 
(NHEK) were purchased from the Clonetics Corp. (San Diego, CA). NHEK cells were grown in 
defined Keratinocyte Growth Medium-KGM™ (Clonetics). 

Cytotoxicity analyses. For toxicity studies, cells were plated at a density of 2,500 cells/cm 2 in 6- 
well plastic tissue culture plates and grown for 24 hr before metal treatment. Triplicate wells of 
cells were treated with increasing concentrations of As, Cr, Cd, Pb, or a mixture of the four 
metals for 24 hr. Concentrations used for As and Cr were 0.3, 1, 3, 10 , and 100 mM. The Cr 
stock was made by mixing a 1:1 ratio of trivalent and hexavalent chromium. Cd and Pb were 
administered at 3, 10, 30, 100, and 300 mM. To prepare the four metal mixture, a 50X stock 
solution was made which, when diluted to its final IX concentration in cell culture, contained the 
levels of each As, Cr, and Cd giving 50% cell killing in individual cytotoxicity assays (LD 50 ). 
The concentrations of these three metals in the IX mixture solution was 7.7 mM As, 4.9 mM Cr, 
6.1 mM Cd. The lead acetate concentration in the IX solution was lOOmM, as we were unable to 
get complete killing at any dose tested. This IX solution was serially diluted at a 1 :3 ratio to get 
0.333, 0.1 1 1, 0.037, 0.0123, 0.004, and 0.0014 dilution groups. Double deionized water was used 
as the vehicle control in all cases. After exposure to the individual metals or metal mixture, cells 
were refed with fresh metal-free KGM medium and incubated for 3 days prior to viability 
analysis by the MTT assay. 

MTT assay. The MTT assay was carried out using a modification of the method of 
Mosmann (1983). MTT was dissolved at 5mg/ml in IX PBS [phosphate-buffered saline]. This 
stock solution was filtered through a 0.2 mM filter and stored at 4°C. Immediately before use, the 
stock solution was diluted to 0.5mg/ml with serum-free KGM to make a working solution. 
Working solution (1 ml) was added to each well of cells after aspiration of medium. Cells were 
incubated at 37°C for 3 hr, after which time, the MTT was removed by aspiration. Cells were 
subsequently lysed by addition of 0.5ml of dimethyl sulfoxide (DMSO). The absorbence at 
550nm of samples as well as a DMSO blank was read on a Microplate Autoreader (Bio-Tek 
Instruments, Inc, Winooski, VT). The absorbance of the DMSO blank was subtracted from all 
values. All absorbance values were standardized to the water vehicle control. 

Additivity model. The nonlinear additivity model selected for fitting the single chemical 
data was based on a Gompertz function. 

= a + Y[exp(-exp(-(p 0 +^x { +$ 2 x 2 +P 3 *3 +M4 )))] + & 0) 
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where 

y is the observed response (%viability), 
x\ is the dose of arsenic (mM), 
xz is the dose of chromium (mM), 
xj is the dose of cadmium (mM), 
x 4 is the dose of lead (mM), 

a is an unknown parameter associated with the minimum response, 
y is an unknown parameter associated with the maximum response 

po is an unknown parameter associated with the intercept on the complementary log log scale, 

pi is an unknown parameter associated with the slope of arsenic, 

P 2 is an unknown parameter associated with the slope of chromium, 

p 3 is an unknown parameter associated with the slope of cadmium, 

p 4 is an unknown parameter associated with the slope of lead, and 

e is an unobserved random error term assumed to have mean 0 and constant variance. 

The Gauss-Newton iterative algorithm was used in PROC NLIN in SAS (version 6.12) to find 

parameter estimates under the assumption of a constant variance across the dose range of all four 

compounds. Upon consideration of the range of sample variances for the single chemical data, a 

weighted least-squares criterion was initially considered where the weight of each observation is 

set to the inverse of the sample variance. For convenience in calculating the prediction intervals 

an unweighted analysis is given here. The predicted responses in the two models were very 

similar. 

To compare the observed response at X m to that predicted under the hypothesis of 
additivity, a prediction interval was used at each of the m=l,. . .,M(=7) mixture points where M is 
the total number of mixture points. Let y Xm be the observed sample mean response at the 
combination X m with sample variance Var( y Xm )• Denote the estimated predicted response 
using the model given in (1) by y Xjn with estimated variance Var(j>x w )• Using the approach 
taken by Gennings et al. (1997), a 100(l-a)% prediction interval for the sample mean response at 
X m under the assumption of additivity is given by 

9x m ± ' (l-a / 2;N-1) yj Var (yX m ) + Var (>'X m ) , (2) 
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where N is the total number of observations. If the observed sample mean, y Xm , is not included 
in the prediction interval, then it is reasonable to conclude that the data do not support the 
assumption of additivity. For decreasing dose response curves, if the observed sample mean is 
less than the lower bound of the prediction interval, then synergism can be claimed at X m , i.e., 
lower viability than expected under additivity. If the observed sample mean exceeds the upper 
bound of the prediction interval, then antagonism can be claimed at X m , i.e., greater viability 
than expected under additivity. It is important to note that these prediction intervals are 
piecewise intervals without correction for multiple testing. 
Overall test of additivity: 

An overall test for additivity can be based on testing the hypothesis that the mean 
response under the hypothesis of additivity is the true mean response. The estimated response 
under the hypothesis of additivity is provided by the model given in (1). The estimated response 
for the true mean is provided by the sample means at each mixture group. Define the M- 
dimensional vector d=E(Y add> x)-E(Y x ), which is estimated as S = y add x -y x . Let 9j be the 
vector of model parameters associated with the model given in (1), i.e., 
Gj =[ayPo Pi P2 P3 P4] '-Let B = [Q[ \i\ ^'"V-mY he the combined vector of model 

[~E 0] , 

parameters and mixture point means. Let Q = , where s S is the pxp variance 



. 0 DJ 

covariance matrix of the model parameters given in qi and 

s 2 D = diag[Var(y(\i )) ...Var(y(x M ))] . Using the multivariate delta method, an approximate 
large sample estimate for the variance of 8 is given by s 2 GWG\ where 

G = I I J . For large samples, 5 [GQG'] -1 81a 1 follows a chi square 

\^ r J J 8;s=\,...M;r=\ p+M 

distribution with M degrees of freedom. Then a Wald-type test for testing if the observed 
mixture points agree with the assumption of additivity simultaneously is given by 

MK ' MSE 

which follows an F distribution with M and N-p degrees of freedom for large samples. 

Once the overall test for additivity is rejected, it is of interest to detect and characterize 
departure from additivity at each of the M mixture points. If M>1 then individual tests for 
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departure from additivity and/or multiple 95% prediction intervals suffer from an inflated 
significance level for the family of tests. Hochberg (1990) proposed a refinement of 
Bonferroni's correction to adjust for multiple tests while maintaining an overall significance 
level. His approach is a step-down sequentially rejective procedure which rejects the hypothesis 
(here, of additivity) associated with the m th statistic and all other hypotheses corresponding to 
smaller p- values when 

p-value( m ) < , for m=l,.. .,M. 

M -m + 1 

Our approach is to use Hochberg's correction on the M single degree of freedom tests of 
(3) when M=l. If the m th p-value associated with a mixture point is less than a/(M-m+l) then 
departure from additivity is claimed at all mixture points with smaller p values. The prediction 
intervals given in (2) are used to characterize the type of interaction detected for the points where 
departure from additivity is detected using Hochberg's correction for multiple testing. 
Comparison of fitted curves: 

When the M mixture points are on a fixed ratio ray, these mixture points can be used to 
estimate a dose response curve along the fixed ratio ray using total dose as the measure of total 
exposure. In this case, an extension to the approach above can be made. Of interest is to 
compare the predicted response along a fixed ratio-mixing ray from two sources. Under the 
hypothesis of additivity, the single chemical data are used to estimate an additivity surface. The 
slice of the fitted surface that corresponds to the fixed ratio ray is denoted f A (x; 9j px i). Suppose 
that enough mixture points are experimentally considered along the fixed ratio ray so that an 
independent fit of the mixture data yields an estimate fM(x; 0 2 qxi)- Further, assume that under 
additivity g(qi) = q2- Thus, the additivity hypothesis is equivalent to testing 1 = g(qi) - q2 = 0. 



s vanance-covanance 



* rs, o 1 , 

Define q* = [qf | q2T an d ^ = I q £ J ' w ^ ere s W* is the variai 
matrix of q*. Define X = g(9j) -0 2 . A large sample approximate variance covariance matrix 

dx 



of X is given by s HW*H\ where H = 
additivity for large samples 



. Then under the hypothesis of 

j=\,...,p+q 
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X\HQ. H'\ X 
MSE 



3 approximately F q , N-(p+q)- 



(4) 



Here, the model given in (1) provides f A (x;qi) and 

f M (x; q 2 ) = a* + g*[exp(-exp(-(b 0 * + b,* x + b„*x 2 )))], (5) 
where x is the total dose along the fixed ratio ray. Then, under additivity, 1 = g(qi) - q 2 = 



= 0. Simultaneous confidence band of f A (x;qi)- f\i(x; q 2 ) along the fixed ratio 



a-a 
* 

Y-Y 

Po-Po* 

i=i 

. PiT 

ray: 

Define d(x;q*) = f A (x;qi)- fM(x; q 2 ) as the difference, as a function of dose, between the 
model for the curve under additivity and using only the mixture data. Then 

Var(d(x; 6* )) *d 2 BQ*B', 
dd(x;Q*) 



where B = 



. Following Miller(1981) and Seber and Wild (1989, pg 

Jj=l,...,p+q 

246), a 100(1 -a)% confidence band for d(x;q*) is approximately for all x 



d(x;q*) e d(x;0)± ft~p + q)F a;p+q<N<p+q) ^Diag(Var(d(x;0)) . 

Noting when zero is not included in this confidence band allows one to determine total dose 
levels that depart from additivity. 



RESULTS 

Summary statistics for the single chemical data and the mixture data are provided in the 
Appendix. The additivity model given in (1) was fit to the single chemical data. A plot of the 
observed and predicted responses is given in Figure 13. The observed and model predicted 
responses for the individual metals are graphically provided in Figure 14. Overall, the fit of the 
single chemical data seems adequate. The estimated model parameters and their associated large 
sample standard errors and p values are provided in Table 1. All of the model parameters were 
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significantly different from zero as the p values are all <0.005. The slope parameters associated 
with each of the four metals are negative and significant, indicating that as the dose of the metal 
increases, the % viability decreases significantly. 

Table 1 : Estimated model parameters for the additivity model given in (1) using NHEK cell line. 



Parameter 


Estimate 


Asymptotic 
Standard Error 


P value 


a 


8.76 


2.58 


O.001 


Y 


107.7 


13.0 


<0.001 


Po 


1.87 


0.658 


0.005 


Pi 


-0.197 


0.058 


O.001 


P2 


-0.259 


0.0669 


<0.001 


P3 


-0.177 


0.0557 


0.001 


P4 


-0.0078 


0.0021 


<0.001 



Figure 15 presents the predicted concentration effect curve under additivity for total 
concentrations of the four metals along the ray associated with the LD 5 o mixing ratio. This curve 
is based on the fit of the single chemical data and the definition of additivity. The asterisks 
indicate the observed sample means at the seven dilution points. From this figure, there is 
evidence of departure from additivity at some of the mixture points. In fact, the general trend 
does not seem to agree with the responses predicted under additivity. 

The observed responses at the seven mixture points are provided in Table 2. This table 
also provides the predicted responses under the hypothesis of additivity. Table 3 presents the 
overall test for departure from additivity given in equation (3). The hypothesis of additivity at all 
seven of the mixture points considered is rejected with a p value <0.001 . Table 3 also includes 
the single degree of freedom tests associated with each mixture point. Using Hochberg's (1990) 
correction for multiple comparisons, the three groups with the smallest p values are significantly 
different from that predicted under additivity as p < a/k for k=5,6,7, when the overall 
significance level is set at 5% (or at 1%). 
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Table 3: Test results for testing the hypothesis that the mean of the mixture group is equal to that 



Statistic 


P value* 


Overall F test ((7,200) df) 




179.1 


O.001 



Individual F tests ((1,200) df) 



Dilution 


Statistic 




0.0014 


6.27 


0.013 


0.004 


1.53 


0.217 


0.0123 


0.62 


0.432 


0.037 


37.2 


O.001** 


0.111 


127.6 


O.001** 


0.333 


2.70 


0.102 


1 


16.8 


O.001** 



Using Hochberg's correction for multiple comparisons, 

these groups are associated with a significant interaction using an overall 5% test. 
Those marked with ** are significant with an overall significance level of 1%. 

The direction of departure from additivity can be assessed through consideration of the 
prediction intervals in Table 2 for mixture points associated with significant departure from 
additivity using Hochberg's correction in Table 3. The 95% and 99% prediction intervals fall 
below the observed response for the lowest concentration of 0.0014 dilution. This indicates that 
there is a greater than expected response at very low concentrations. Perhaps this is due to an 
hormesis (i.e., growth stimulating) effect. Further, for the dilutions of 0.037 and 0. 1 1 1 , the 
observed sample means fall below the 95% prediction intervals. That is, the observed responses 
show more cytotoxicity than predicted under additivity. This may be associated with an overall 
conclusion of a cytotoxic synergism at these concentrations. By comparison, the mixture of the 
four LD50s (i.e., dilution is 1.0) show evidence of a less than additivity association in terms of 
cytotoxicity, as the observed mean response was less cytotoxic than predicted under additivity. 
This may be associated with an overall comclusion of a cytotoxic antagonism at this mixture. 

With seven mixture points experimentally evaluated there are enough data to actually fit a 
concentration effect curve to these data and compare the fitted function to that predicted under 
additivity. Preliminary evaluation of the mixture data indicated that the relationship was more 
quadratic than linear on the complementary log-log scale. For this reason, the model given in (5) 
was used to fit the mixture data which includes a quadratic term in concentration. Figure 15 
presents the fitted mixture data concentration effect curve (dashed line) plotted with the 
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predicted response curve under the hypothesis of additivity (solid line) and the observed mixture 
points (asterisks). The overall test for a difference between the fitted curve from the mixture 
data and the curve predicted under additivity as given in (4) was rejected (p<0.001). 

From Figure 15, there seems to be an increasing departure of the two curves in the 
concentration range of roughly 10 to 40 total concentration units. This is depicted in Figure 16 
as the difference between the predicted response of the mixture data and that predicted under 
additivity. The 95% simultaneous confidence band on the difference between the curves (Figure 
16) does not include zero in the range from 8 to 36 mM of the mixture. In this range, the 
difference is significantly negative indicating the response of the mixture was associated with an 
increase in cytotoxicity as compared to that predicted under additivity, i.e., synergism. In 
addition, the difference is positive in the range from 80 to 120 mM indicating a decrease in 
cytotoxicity as compared to that predicted under additivity, i.e., antagonism. Thus, from this 
example, we have evidence of different characterizations of interaction along the same fixed 
ratio mixture of the four metals. 
DISCUSSION 

The overall seven degree of freedom test for departure from additivity indicated that for 
at least one of the seven mixture groups there is a difference between the mean response as 
estimated by the additivity model and the unrestricted mean. The former is estimated by the 
model given in (1) and the latter by the sample means. Similar to the post hoc tests used in 
analysis of variance, Hochberg's step-down sequentially rejective procedure results in at least 
three concentration groups being associated with a significant difference from additivity. 

The fit of the additivity model assumed a common estimate for a, g, and bo. This restricts 
the maximum and minimum values for the range of the response to be the same: a minimum 
value of percent viability of 8.8% and a maximum value of 101 .3%. The lowest mixture 
combination had a mean response of 1 16.6. This is clearly above that observed in the single 
chemical data. It is likely that this enhancement of cell viability at the lowest mixture level is 
indicative of the presence of "hormesis." Hormesis, defined by Stebbing (1982) as the 
stimulatory effects caused by low levels of toxic agents, was originally developed under different 
terminologies over a century ago. Its origin, scientific development, historical perspectives, and 
modern implications have been ably reviewed by a number of contemporary scientists (Stebbing, 
1982; Calabrese, 1997; Calabrese and Baldwin, 1997a,b; Stebbing, 1997). Since we are 
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observing growth stimulatory effects at the lowest level of the metal mixture, this might be the 
first observation that a mixture of As, Cd, Cr, and Pb caused a hormesis effect. The interactive 
effect changes from additive cytotoxicity to synergistic cytotoxicity in the range of about 8 to 36 
mM of the metal mixture (Figure 4). It should be noted that when we talk about synergism or 
antagonism here we are referring to cytotoxicity; it is completely opposite if we refer to cell 
viability. Thus, we observe hormesis at the lowest mixture level and synergistic cytotoxicity at 
higher levels (roughly 8 to 36 mM). What we are seeing is a typical b-curve (or J- or U-shaped 
curve) for the growth hormesis (Calabrese, 1997). As Stebbing (1997) proposed, hormesis may 
be the cumulative consequence of transient and sustained over-corrections by feedback growth 
regulatory mechanisms to low levels of inhibitory challenge. When the concentration level 
increases, toxicity becomes more and more prevalent and the growth regulatory mechanisms are 
overwhelmed. Growth recovery to control level is thus impossible. At this stage, a number of 
mechanisms of toxicity for the different metals may work in concert, antagonism toward cell 
viability (or synergistic cytotoxicity) is therefore observed (Figures 15 and 16). Interestingly, at 
the highest mixture concentration tested in NHEK (119 mM), we again observed antagonistic 
cytotoxicity. We have also seen this same trend in immortalized keratinocyte cell lines when 
treated with a similar metal mixture (Bae et al., 2001). One hypothesis for why this is occurring 
is that at very high metal concentrations, cell protective mechanisms are induced. We are 
currently exploring this possibility by quantitation of the molecules glutathione and 
metallothionein; both of which have been shown to confer protection against metal-mediated 
cytotoxicity. 
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APPENDIX 



NHEK Cell line 



Chemical 


Dose 


Sample Mean 


Sample 


Sample Size 








Variance 




Arsenic 


0 


100.0 


45.7 


9 




0.3 


108.9 


146.6 


9 




1.0 


99.6 


78.1 


9 




3.0 


84.1 


279 


9 




10.0 


45.1 


268 


9 




30.0 


8.4 


108 


9 


Chromium 


0 


100.0 


137 


6 




0.3 


105.1 


119 


9 




1.0 


99.7 


124 


9 




3.0 


79.8 


304 


9 




10.0 


23.9 


446 


9 




30.0 


6.6 


78 


9 


Cadmium 


0 


100.0 


27 


6 




3 


92.6 


.130 


9 




10 


51.9 


1236 


9 




30 


20.2 


829 


9 




100 


3.8 


23 


9 




300 


4.6 


28 


9 


Lead 


Q 


100 0 


8 4 






3 


95.5 


79 


9 




10 


91.0 


63 


9 




30 


98.8 


236 


9 




100 


93.7 


1838 


9 




300 


28.9 


297 


9 


LD50 Mixture * 


Total 








(dilution) 


Concentration 








0.0014 


0.16 


116.6 


315 


9 


0.004 


0.48 


95.7 


1.9 


9 


0.0123 


1.46 


96.0 


216 


9 


0.037 


4.39 


76.5 


92 


9 


0.11 


13.2 


64.3 


14 


9 


0.33 


39.6 


50.1 


262 


9 


1.0 


118.7 


31.1 


207 


9 



* The predetermined LD 50 s were 7.7, 4.9, 6.1, 100 mM for arsenic, chromium, cadmium, and 



lead, respectively. 
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Example 6. Experimental Designs for Mixtures of Chemical s Along Fixed Ratio Rays 
Abbreviations: 



As: Arsenic 

ATSDR: Agency for Toxic Substances and Disease Registry 

5 Cd: Cadmium 

CHEM: Chemical 

Cr: Chromium 

ED 50 : Dose/concentration associated with 50% response 

GLM: General linear model procedure in SAS 

10 LD 50 : Dose/concentration associated with 50% lethality 

Pb: Lead 

RSM: Response Surface Methodology 



Introduction 

Determining and characterizing the nature of interactions among components of a 

1 5 combination of c drugs or chemicals is a problem of current interest (where c is the number of 
drugs/chemicals in a mixture). Although assessments based on single drug/chemical exposure 
enable us to acquire fundamental knowledge about individual drugs or chemicals under carefully 
controlled conditions, they do not reflect real-world exposures. Thus, it is often of interest to 
study the effects of exposure to multiple drugs/chemicals. Of ultimate interest in such studies is 

20 the determination and characterization of interactions among the components in a mixture. For 
example, Gennings et al. (1) report on a study of the nature of the interaction involving the 
mixture of four metals. The four metals chosen for the study were arsenic (As), cadmium (Cd), 
chromium (Cr), and lead (Pb), which are among the top contaminants in site frequency count by 
the Agency for Toxic Substances and Disease Registry (ATSDR) Completed Exposure Pathway 

25 Site Count Report (2). In addition, human health risk assessment associated with exposure to 
disinfection by-products (DBPs) in drinking water is of concern because of the wide spread 
exposure of persons who receive disinfected water. Other examples of human exposure to 
combinations of agents can be found in the treatment of numerous diseases including cancer, 
AIDS, diabetes and asthma. These examples illustrate the importance of studying 

30 mixtures/combinations of drugs or chemicals. Determining departures from additivity for a 

combination of drugs or chemicals is a problem that has been considered by many authors (3-7). 
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Classical Methodology for Detecting and Characterizing Departures from Additivity 
2.1 Isobolograms 

The classical method for detecting and characterizing departures from additivity between 
combinations of drugs or chemicals is the isobologram. The isobologram, introduced as a 
5 graphical tool by Fraser (8,9), is a plot of a contour of constant response of the dose-response 
surface associated with the combination superimposed on a plot of the same contour under the 
assumption of additivity. Their use was extended by Loewe and Muischnek (10), Loewe (1 1), 
and Berenbaum (12) and reviewed by Gessner (13), Wessinger (14), and Berenbaum (15). For a 
two-component mixture, the analysis of an isobologram compares the observed isobol (e.g., 

10 combination ED 50 ) to the line of additivity. The line of additivity is formed by joining the ED 50 
associated with each of the individual components calculated from the dose response data for the 
individual components. Figure 17 presents illustrations of possible isobolograms for a 
combination of two drugs/chemicals. As indicated, if the isobol is below the line of additivity, a 
synergism is claimed. On the other hand, if the isobol is above the line of additivity, an 

15 antagonism is claimed. However, there are shortcomings associated with the use of 

isobolograms. For instance, the method used in the construction of an isobologram typically 
does not take data variability into account. Additionally, since it is a graphical method, 
isobolograms effectively are limited to the study of combinations of two or three drugs or 
chemicals. 

20 2.2 Interaction Index 

The interaction index, introduced by Berenbaum (12), provides a convenient method to 
determine and characterize departures from additivity for a combination of c> 2 or 3 
components. The interaction index, II, is defined by 

25 // = ^ + £ + ... + ^ (1) 

ED mfi {CHEM x ) ED l00fX (CHEM 2 ) ED l00fl (CHEM c ) 

where c is the number of components, X\, X2,. . .JC C are the doses in combination associated with 
a desired effect, and EDwo^CHEMj), i=\,...,c is the dose of the rth component that, when 
administered alone, produces the same effect. When the interaction index, defined in equation 
(1 ), is equal to 1 the c components interact additively; when II is greater than 1 the components 
30 interact antagonistically; and when II is less than 1 the components interact synergistically. 
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Again, it should be noted that the individual component dose-response information is required to 
calculate the interaction index. As described by Berenbaum (12), the interaction index is directly 
related to the isobologram, i.e., when n=l, the isobol is coincident with the line of additivity; 
when II > 1, the isobol bows above the line of additivity; and when II < 1, the isobol bows below 
5 the line of additivity. An advantage of using the interaction index over the isobologram is that 
the interaction index is not limited to combinations/mixtures of just two or three components. 
However, as developed by Berenbaum (12), the biological variability associated with the data is 
not taken into account by the interaction index. 
2.3 Statistical Models 

10 

Statisticians frequently use models of the form 

c c c c c c 

y = fa + 2 Pi* + 2 2 fy x t x j + 2 22 - + 012.^2-^ 

,=l ,=1 y = l ,=i y=iifc=i 

i<j i<j<k 

15 to approximate the relationship between a response of interest, Y, and concentrations of c 
chemicals (x\j2,- • -yX c )- 

Carter et al. (4) showed that a relationship exists between the interaction index proposed 
by Berenbaum (12) and the parameter in a statistical model that is associated with the interaction 
of the components of the combination. Without loss of generality, consider that the 

20 combination/mixture of interest involves two chemicals and that the response is continuous. 
Therefore, following the logic of Carter et al. (4), for the linear models case, the relationship 
between the response and the doses or concentrations of the components in combination can be 
expressed as 

ft = 00 + A*l + 02*2 + #2*12 > (2) 

25 where 

ju = is the mean response, E(Y), 
J3q is the unknown intercept, 

J3\ is the unknown slope parameter associated with the first component, 
pi is the unknown slope parameter associated with the second component, 
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/?i2 is the unknown parameter associated with the interaction of the two components, 
and x\ and xi are the doses of the respective chemicals. 

From the model defined in equation (2), the EDm^CHEM,) for the respective 
components can be derived to be 

M~ fio 



ED mtx {CHEM x ) = 



fix 
M-fio 



ED m „(CHEM 2 ) 

Pi 

Thus, after algebraic manipulation, the model defined in equation (2) becomes 

01*1 + 02*2 + #2*12 = ! 
M-fio M-fio M-fio 



i _ A 2*1 2 _ *1 , x 2 

»-p Q (M-fio)/ (M-fio)/ ' 
/ fii /Pi 

10 Therefore, it follows that when /?i2=0, the combination of components 1 and 2 is additive, i.e., 
the isobologram is coincident with the line of additivity and the interaction index equals 1. 
Similarly, when /?i 2 >0, the combination of components 1 and 2 is synergistic, i.e., the 
isobologram bows below the line of additivity and the interaction index is less than 1; and when 
fiuO, an antagonism is present, i.e., the isobologram bows above the line of additivity and the 

15 interaction index is greater than 1. This demonstrates the algebraic equivalence between the 
statistical model and the interaction index. Gennings et al (16) demonstrated the experimental 
convergence of the statistical modeling approach and the interaction index. The number of 
components that can be considered in the statistical model can be generalized to c, and data 
variability is appropriately accounted for in the resulting inference. 

20 The various methods used to determine and classify departures from additivity described 

above utilize both single-compound data as well as combination data. Consider the situation in 
which the single-compound dose-response data are not available. The approach discussed in the 
following sections of this paper allows one to test the null hypothesis of additivity using only 
combination/mixture data collected along a fixed ratio ray. Thus, single-compound dose- 

25 response data are not needed. 
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New Methodology 

Problems with statistical modeling are associated with the size of the experiment required 
to generate data to support the model. Factorial experiments, e.g. 2 C or 3 C are often considered. 
When c is large such experiments may not be feasible, so, in a sense, this approach is limited to 
5 combinations of relatively few drugs or chemicals. An alternative to the traditional factorial 
design for studying interactions, and the design to be considered here, is the ray design. Ray 
designs, described by Martin (17), Mantel (18), Finney (19), Bruden and Vidmar (20), and 
others, are used to study mixtures of c drugs or chemicals at a fixed mixing ratio, [aj :a 2 : . . . :a c ], 
c 

where £ a t = 1 , with the total "dose", t, varying. The fraction of the total "dose" represented 
/=i 

10 by the rth drug or chemical is a, and the amount of the rth drug or chemical in the mixture is a,*. 
This approach is appealing since the dimensionality of the study is reduced along each ray, i.e., 
each ray can be considered as an individual drug or chemical with only the total "dose" varying. 
For example, in a study involving c drugs or chemicals the fitted model based on a response 
surface approach is a (c+l)-dimensional surface. In contrast, the fitted model based on a ray 

15 design defines a set of 2-dimensional dose response curves. 

What can be stated about departures from additivity in the mixture? Meadows (21) 
showed that when mixture data are collected along a fixed ratio ray, the additivity model reduces 
to a simple linear regression model. In addition, the interaction model reduces to a higher order 
polynomial model. Thus, the test for additivity is equivalent to the test of the adequacy of the 

20 simple linear regression model. Consider that the combination/mixture of interest involves c 
drugs or chemicals and that the response of interest is continuous. The underlying additivity 
model, i.e., the model with no cross-product terms is defined by 

r = + + &* c > (3) 

25 

where 

Y is the observed response, 
Xi is the dose of the rth drug or chemical, 
fto is an unknown parameter associated with the intercept, and 
30 Pi is an unknown parameter associated with the slope of the rth drug or chemical. 
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When the mixing ratios are invoked, the dose of the z'th drug or chemical is x, = a,/, where 
a, is the mixture fraction for the z'th drug or chemical and t is the total dose. As a result, the 
additivity model becomes 
5 Y = /? 0 + fat + Pi^it + - + Pc"ct 

= Po + (M- + £2*2 + - + P c a c )t 

= Po + Pit, (4) 

where /?,* = ^ /^a,- . For convenience, we assume the experimental region along the fixed 
1=1 

ratio ray in terms of total dose is transformed to the region -1 < t < 1. Thus, under additivity, the 
dose response relationship along the ray can be described with a simple linear regression on total 
10 dose. 

It also follows that when the slope of the regression line for total dose is 
M + Pz a 2 + - + Pc a c> 
the interaction index, defined in equation (1) equals 1. This would suggest that the single- 
component dose response data would be required. In the absence of single chemical data, the 

c 

15 slope bj for each drug/chemical alone is unknown so that the hypothesis H 0 : fi\ = ]T fifa can 
not be directly tested. However, consider the following model 



c c c c c c 

Y = Po + YuPi X i + HHPij X i X j + E E HPijk x i x J X k+ - + P\2... c X\X2-Xc 

, = 1 , = 1 j = \ ,=1 jmlkml 

i<j i<j<k 

Notice that this model is the model that would be supported by a factorial experiment and the 
20 PipPijh ■ ■ -Pi2...c terms are coefficients associated with the various two, three-factor and higher 
order interactions. The ray design will not support this model; however, invoking the mixing 
ratio associated with the ray design results in 
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V( = l J ; = I y = i 

V 1< J ) 



c c c 




+ ... + (^...a^ .Jt c 



= p 0 + p[t + p\t 2 + p\t l + ... + p*t c 



(5) 



It follows that interactions among pairs of chemical components are associated with second- 
5 degree terms, interactions among three chemicals are associated with third degree terms, etc. 

Of ultimate interest is the determination of departure from additivity among a particular 
combination/mixture of drugs or chemicals. When comparing the model under additivity to the 
interaction model, defined in equations (4) and (5), respectively, evidence of curvature indicates 
departure from additivity; i.e., there is interaction among the compounds if at least one P , ^ 0, 
10 i'=l , . . . ,c. Thus, any polynomial lack of fit associated with the additive model 



would be associated with a lack of additivity. Meadows (21) showed that the test statistic for the 



Y = Po + PU 



null hypothesis of additivity, H 0 : 



Pi 



= 0 , is given by 



X. 



(6) 



15 



4. 



Experimental Design 

Experimental design implications for studying a c component mixture include the 



following: 



20 



(1) Place a minimum of c+1 points on the ray of interest so as to maximize the power of 
the test for lack of fit of the additivity model, i.e., Y = p 0 + p x t . 

(2) Replicate the experiment at these points to make the lack of fit test possible. 
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When the response variable is continuous and the method of least squares has been used 
to estimate the model parameters, Meadows (21) showed that we can incorporate the statistical 
results of Jones and Mitchell (22) to determine values of total dose that maximize the design's 
5 ability to detect lack of fit or departure from additivity. 

The overall lack of fit answers the question as to whether or not there is a departure from 
additivity. Rejection of his hypothesis that simultaneously tests that the interaction parameters 
are equal to zero, i.e., Ho: Pi=fiy=. • = A = 0, implies that interaction is present among the 
chemicals globally. Thus, if the overall test for additivity is rejected, tests of the form 
10 H 0 : fij = 0 

Hi:/?* * 0 

using Hochberg's (23) correction for multiple testing, can be used to answer the question of 
whether or not a y-factor interaction exists. If such an interaction is detected recall that 

ti = tt WjPy A = I £ t *flj*kh» etc - 

i=l 7=1 i=l j=lk=\ 

i<j 

1 5 Here interest will be focused on which of the y-factor interactions are present. This can be 

determined by performing ^ j additional ray experiments at the same ratio as were present on 
the original ray. 

5. Illustration 

20 The methodology introduced in this paper is illustrated with cytotoxicity data obtained 

from assessing interactions among arsenic (As), cadmium (Cd), chromium (Cr), and lead (Pb) in 
human keratinocytes. The experimental data were obtained from Ray Yang and colleagues at 
Colorado State University. The endpoint of interest is the percent viability of treated NHEK 
cells using the MTT assay. The mixture point of interest for As, Cr, Cd, and Pb contained the 

25 estimated LD 50 s of 7.7uM, 4.9 uM, 6.1 uM, and 100 uM, respectively. This IX solution was 
serially diluted at a 1:3 ratio to get 0.333, 0.1 1 1, 0.037, 0.0123, 0.004, and 0.0014 dilution 
groups. Double deionized water was used as the vehicle control in all cases. 
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After exposure to individual metals or metal exposure, cells were refed with fresh metal- 
free KGM medium and incubated for 3 days prior to viability analysis by the MTT assay. 
Details of the experimental protocol and methods are described elsewhere (24, 1) and are not 
included here. The summary statistics for the LD 50 mixture data presented in Table 1 are 
5 linearized cytotoxicity response data from Gennings et al. (1). 



Table 1. Summary Statistics for the LD50 Mixture Data Based on the Linearized Cytotoxicity 



Mixture 


Total Dose 3 


Sample Mean 


Sample 


Sample size 


Dilution 


(uM) 




Variance 




0 


0 


1.81 


0.23 


9 


0.0014 


0.2 


2.24 


0.20 


6 b 


0.004 


0.5 


1.65 


0.55 


9 


0.0123 


1.5 


1.58 


1.07 


8 b 


0.037 


4.4 


0.77 


0.11 


9 


0.111 


13.2 


0.40 


0.01 


9 


0.333 


39.6 


0.02 


0.18 


9 


1 


118.7 


-0.55 


0.27 


9 



a LD 50 mixing ratio (7.7pM, 4.9uM, 6.1|iM, and lOOuM) for As, Cr, Cd and Pb, respectively. 
1 0 b Overall, four data values are missing due to the transformation on the response. 



The nonlinear additivity model selected for fitting the single compound data by Gennings 
et al. (1) was based on a Gompertz function where the mean viability was modeled as 

ju = a + y exp(-exp(-(/? 0 + ^ /?,*, ))) . From this model, a is the parameter associated with the 

i=i 

1 5 minimum mean response and g is the range of mean response values. Therefore, for this 

example, it is reasonable to assume that a transformation on the response, conditioning on the 
values of a=8.76 and y=109 obtained by Gennings et al. (1), will induce linearity in the 
additivity model. As a result, the additivity model becomes 

-Sod -lo/^ 



(7) 



20 As shown in Section 3, since the mixture data were collected along a fixed ratio ray, the 
additivity model can be re-written as 



Pit- 



Mm.add = - lo S ~ l0 § 
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Additionally, the interaction model along the same fixed ratio ray becomes 

^.interaction = ~ logf " logf^-^-N = 0 O + fit + /? 2 V + + /? 4 V . (8) 



7 

Therefore, conditioning on the values of cc=8.76 and y=109, the transformation 

on the observed responses for the mixture data was performed. The 



log| - log^- 



y - 8.76 
109 



additivity model given in equation (7) and the interaction model given in equation (8) was fit to 
the mixture data using the method of least squares. The GLM procedure of SAS® was used to 
estimate the unknown parameters in equations (7) and (8). Parameter estimates and their p 
values are provided in Table 2. Figure 18 presents the fitted concentration effect curve under 

10 additivity for total concentrations of the four metals along the ray associated with the LD 5 o 

mixing ratio. The asterisks (*) indicate the observed transformed responses at the seven dilution 
points. From this figure, there is some question as to whether the data fall along the line of 
additivity. In comparison, Figure 19 presents the observed mixture data and the fitted interaction 
(higher order polynomial) model. The dots (•) indicate the design locations of the total dose 

1 5 values selected by the A i -optimal design, which are presented by total doses of 0, 1 6, 59.4, 

102.7, and 1 18.7 uM. Notice that the values selected as the Ai-optimal design are symmetrically 
spread out throughout the total dose region, whereas the majority of the points used in the current 
study are directed toward the lower total dose region. An enlarged version of the lower total 
dose region of the plot of the fitted interaction model is presented in Figure 20. 

20 

Table 2: Estimated Model Parameters for the Additivity Model Given in Equation (7) and the 



Parameter 


Estimate 


S.E. 


P value 


Additivity Model: 








Po 


1.36 


.12 


< .001 




-0.02 


2.34 x 10" 3 


<.001 


SSres=38.04, dfRE S =66 




Interaction Model: 








fa 


1.96 


.12 


<001 


fit® 


-0.36 


.09 


<001 


#(t 2 ) 
#(t 3 ) 


0.03 


.01 


.010 


-0.0006 


2.44 x 10" 4 


.024 
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Mt 4 ) 



0.00000319 1.43 x 10 -6 
SSres=20.49, dfREs=63 



.030 



The test statistic for the null hypothesis of additivity, H 0 : 



0, is given by 



_ _ R(o 2 | e x )i(c - 1) _ [ss Rcs> 

reduced ~ ^Res,full)/ 3 

F = 7 = 7 

_ (38.04 - 20.49)/3 
0.3252 

= 17.99 

5 Table 3 presents the overall test for departure from additivity given in equation (6). 



Table 3: Test results for testing the hypothesis of additivity, as well as the hypotheses that the 



Hypothesis 


F 


P value 


Overall test for Additivity: 




H 0 : 0f=/h=fr=O 


17.99 (3, 63) 


<0.001 


Individual tests: 






H 0A : A=0 
H ob : #=0 
Hoc: A=0 


23.95 (1,65) 
16.54(1,64) 
4.92(1,63) 


O.OOl** 
O.001** 
0.0302* 



* Using Hochberg's correction for multiple comparisons, these tests are associated with a 
significant ./'-factor interaction using an overall 5% test. Those marked with ** are significant 
with an overall significance level of 1%. 



i on this test, we reject the null hypothesis of additivity (p-value < 0.001) and conclude that 
at least one of the y-factor interactions exists, j=2,. . .4. Since the overall test for additivity is 
1 5 rejected, it is of interest to determine whether two, three, or four-factor interactions are present. 
Therefore, we want to test the following hypotheses using Hochberg's (23) correction for 
multiple testing: 
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H 0A :#=° vs. H 1A :#*0 

H 0B :/? 3 '=0 vs. H IB :#*0 

H oc : /?* = 0 vs. H 1C : p\ * 0 
Table 3 also presents the single parameter tests associated with each of the y'-factor interactions, 
j=2,. . .4. Using Hochberg's (23) correction for multiple comparisons, all three parameters are 
significantly different from zero, when the overall significance level is set at 5%. However, if 
5 we consider the case where an overall significance level of 1% is used, only the two smallest p- 
values are significant using Hochberg's correction. Therefore, we conclude that a 3-factor 
interaction exists implying that the 2-factor interactions are not constant. Now it is of interest to 
determine which three metals are interacting with one another. This can be accomplished by 

performing = 4 additional experiments at the same ratios of metals that were used in the 

10 original ray. Table 4 gives the ratios of compounds, along with the corresponding total dose, to 
be used for the four additional experiments. This approach limits the inferences of the original 
experiment, as well as the four additional experiments, to be made about the particular mixing 
ratio used in the experiment. 



15 Table 4: Ratios of compounds to be used for the four additional experiments, which are based on 



the LD50 mixing ratio (TJuM, 4.9uM, 6 luM, and lOOuM). 





As 


Cr 


Cd 


Pb 


Total dose 


4-factor combination: 

Original Experiment 


6.5% 


4.1% 


5.1% 


84.3% 


118.7 uM 


3-factor combination: 
Experiment #1 


41.2% 


26.2% 


32.6% 




18.7 uM 


Experiment #2 


6.8% 


4.4% 




88.8% 


112.6 uM 


Experiment #3 


6.8% 




5.3% 


87.9% 


113.8 uM 


Experiment #4 




4.4% 


5.5% 


90.1% 


lll.OuM 
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Conclusion 

It was shown that the classical methodology used in evaluating an interaction 
requires single drug/chemical data. In Section 3 it was shown that the evaluation of 
interactions could be accomplished with a ray design that did not generate single 
drug/chemical data. When a ray design is used, departure from additivity is associated 
with higher order polynomial terms in a linear model. Additivity, or absence of 
interaction, is described by a simple linear model in terms of total dose. As a result, we 
have shown that we can obtain information about departures from additivity from data 
collected along a fixed ratio ray. This result is important in that it permits a reduction in 
the total experimental effort for studying a combination when compared to that associated 
with a traditional factorial design. Additionally, by incorporating the approach taken by 
Jones and Mitchell (22), we have presented methodology for determining optimal levels 
along the fixed ratio ray (total dose) to be considered in the experiment for detecting 
model inadequacy. 



REFERENCES 

1 . Gennings C, Carter WH Jr., Campain J, Bae D-S, Yang RSH. Statistical Analysis 
of Interactive Cytotoxicity in Human Epidermal Keratinocytes Following Exposure to 
a Mixture of Four Metals. Journal of Agricultural, Biological, and Environmental 
Statistics 7(l):58-73 (2002). 

2 . ATSDR. 1997 CERCL A Priority List of Hazardous Substances That Will be the 
Subjects of Toxicological Profiles & Support Document. Agency for Toxic 
Substances and Disease Registry, U.S. Department of Health and Human Services, 
November, 1997. 

3 . Gessner PK and Cabana BE. A Study of the Interaction of the Hypnotic Effects 
and of the Toxic Effects of Chloral Hydrate and Ethanol. Journal of Pharmacology 
and Experimental Therapeutics 174: 247-259 (1970). 

4. Carter WH Jr, Gennings C, Staniswalis J, Campbell E, and White K. A Statistical 
Approach to the Construction and Analysis of Isobolograms. Journal of the American 
College of Toxicology 7: 963-973 (1988). 



142 



5. Kelly C and Rice J. Monotone Smoothing With Application to Dose-Response 
Curves and the Assessment of Synergism. Biometrics 46: 1071-1085 (1990). 

6. Gennings C, Carter WH Jr. Utilizing Concentration-Response Data From 
Individual Components to Detect Statistically Significant Departures From Additivity 
in Chemical Mixtures. Biometrics 51:1264-1277 (1995). 

7. Dawson KS, Carter WH Jr., and Gennings C. A Statistical Test for Detecting and 
Characterizing Departures from Additivity in Drug/Chemical Combinations. Journal 
of the Agricultural, Biological, and Environmental Statistics 5: 342-359 (2000). 

8. Fraser TR. An Experimental Research on the Antagonism Between the Actions of 
Physostigma and Atropia. Proc. R. Soc. Edinburg 7: 506-511 (1870-1871). 

9. Fraser TR. The Antagonism Between the Actions of Active Substances. British 
Medical Journal ii: 485-487 (1872). 

10. Loewe S and Muischnek H. Uber Kombinations-Wirkungen. 1. Mitteilung: 
Hilfsmittle der Fragestel-lung. Archives of Experimental Pathology and 
Pharmacology 114: 313-326 (1926). 

1 1 . Loewe S. The Problem of Synergism and Antagonism of Combined Drugs. 
Arzneimittle forshung 3: 285-290 (1953). 

12. Berenbaum MC. Criteria for Analyzing Interactions Between Biologically Active 
Agents. Journal of the American Statistical Association 78: 90-98 (1981). 

13. Gessner PK. The Isobolographic Method Applied to Drug Interactions. In: Drug 
Interactions (Morselli PL, Garattini S and Cohen S, ed). New York: Raven Press, 
1974; 349-362. 

14. Wessinger WD. Approaches to the Study of Drug Interaction in Behavioral 
Pharmacology. Neurosci. Behav. Rev. 10: 103-113 (1976). 

15. Berenbaum MC. What is Synergy? Pharmacol. Rev. 41 : 93-141 (1989). 



16. Gennings, C, Carter, W.H., Jr., Campbell, E.D., Staniswalis, J.G., Martin, T.J., 
Martin, B.R., White, K.L. Jr.: Isobolographic Characterization of Drug Interactions 
Incorporating Biological Variability. Journal Pharmacology and Experimental 



143 



Therapeutics , 252:208-217 (1990). 



17. Martin JT. The Problem of the Evaluation of Rotenone-Containing Plants. VI. 
The Toxicity of 1 -elliptone and of Poisons Applied Jointly, with Further 
Observations on the Retonone Equivalent Method of Assessing the Toxicity of 
Derris Root. Annals of Applied Biology 29: 69-81 (1942). 

18. Mantel N. An Experimental Design in Combination Chemotherapy. Annals of 
the New York Academy of Sciences 76: 909 (1958). 

1 9. Finney DJ. Probit Analysis. London: Cambridge University Press, 1 964, 1971. 



20. Bruden M and Vidmar T. Optimal 3x3 Factorial and 3-ray Designs for Drug 
Interaction Experiments with Dichotomous Responses. Communications in Statistics 
18(8): 2835-2859 (1989). 

21 . Meadows SL. Design Optimality And Sample Size Planning For Mixtures Studies 
[PhD Thesis]. Richmond, VA: Virginia Commonwealth University, 2001. 

22. Jones ER and Mitchell TJ. Design Criteria for Detecting Model Inadequacy. 
Biometrika 65: 541-551 (1978). 

23. Hochberg Y. A Sharper Bonferroni Procedure for Multiple Tests of Significance. 
Biometrika 75: 800-802 (1988). 

24. Yang RSH, ed. Toxicology of Chemical Mixtures: Case Studies, Mechanisms, 
and Novel Approaches. San Diego: Academic Press, 1999. 



144 



Example 7. Analysis of Mixtures of Drags/Chemicals Along a Fixed Ratio Ray Without 
Single Chemical Data to Support an Additivity Model 
1. Introduction 

The determination and characterization of departures from additivity among drugs 
5 or chemicals in mixtures/combinations is often of importance. The definition of 

additivity assumed throughout is given by Berenbaum (e.g., 1985) and is equivalent to 
that based on the classical isobologram for the combination of two chemicals (e.g., 
Loewe and Muischnek, 1926; Loewe, 1953). That is, in a combination of c chemicals, let 
Ej represent the concentration/dose of the i th component alone that yields a fixed 
10 response, mo, and let xj represent the concentration/dose of the i th component in 

combination with the c agents that yields the same response. According to this definition 
of additivity, if the substances combine with zero interaction, then 

which is the interaction index considered by Dawson et al. (2000). A response surface 
15 that satisfies this definition of additivity is termed an 'additivity surface'. The additivity 
surface can be used to test the null hypothesis that the chemicals combine in an additive 
fashion. The experimental data necessary and sufficient to support the estimation of the 
additivity surface are single chemical dose-response data (e.g., Gennings, et al., 1997; 
2000). 

20 Also of importance is the consideration of the experimental design used when 

studying mixtures/combinations of drugs or chemicals. Factorial experiments involving c 
drugs/chemicals, such as 2° or 3°, are often considered. However, when c is large, such 

experiments may not be feasible. Thus, use of such factorial designs is limited to 
combinations of relatively few drugs or chemicals. Alternatively, ray designs, described 
25 by Martin (1942), Mantel (1958), Finney (1964), Bruden and Vidmar (1989), and others, 
are used to study mixtures of c drugs or chemicals at several levels of a fixed mixing ratio 

c 

[ar. a 2 :...:a c ], such that £ a t = 1. 

i=i 

The proportion of the total "dose" represented by the z'th chemical (i=l,2,...,c) is 
denoted by a, and the amount of the ith drug or chemical in the mixture is ait, where t is 
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the total dose. The mixture can then be experimentally observed at several dilutions, 
creating a fixed-ratio ray. Ray designs can be used to support the estimation of a 
response surface (e.g., Bruden and Vidmar, 1989). 

Alternatively, inference may be focused to a particular relevant fixed-ratio ray. A 

5 mixing ratio may be considered relevant, for example, if the mixture under study is 
observed environmentally. This approach is appealing since the dimensionality of the 
study is reduced along the ray of interest. The drug/chemical combination along the ray 
can be considered as an individual drug/chemical with only the total "dose" varying. For 
example, in a study involving c drugs or chemicals the dose-response relationship is 

10 described by a (c+l)-dimensional surface. In contrast, the dose-response relationship on 
a given fixed ratio ray can be described by a 2-dimensional dose-response curve in the 
total dose, t. The methodology developed in this paper assumes that the mixture data are 
collected along a fixed ratio ray. 

In the data considered here, there are three mixture points along a fixed ratio ray. 

1 5 Following the approach described in Gennings et al. (2002), a generalized linear model is 
fitted to the mixture data and then the estimated dose-response curve (over the total dose 
range) is statistically compared to a dose-response curve fit using the additivity model. If 
the two curves are significantly different, then departure from additivity can be claimed 
along the fixed-ratio ray. Otherwise, evidence of departure from additivity cannot be 

20 claimed, providing evidence that additivity is a reasonable assumption for risk 

assessment. This general approach is described in detail by many authors, including 
Berenbaum (1985), Kelly and Rice (1990), Gennings and Carter (1995), Gennings et al. 
(1997; 2000). 

However, for the study considered herein, comparing the positive control values 
25 from the mixture studies to that predicted at the same dose levels from the original single 
chemical studies indicated that the dose-response curves shifted. Thus, the single 
chemical data and hence the additivity surface cannot be used to compare to the mixture 
data. Differences between the two curves may be due to changes in the experimental 
environment and not due to a true interaction among the test chemicals. The purpose of 
30 this paper is to develop an approach that can be used to test for interaction among c 
chemicals when single chemical data are either unreliable or are not available. 
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The mixture selected for the study used to illustrate the methodology proposed 
here is that of four trihalomethanes formed as disinfection by-products (DBPs) when 
water is chlorinated. The four components are bromodichloromethane (BDCM), 
chlorodibromomethane (CDBM), chloroform (CHCU), and bromoform (CHBr3). The 
5 endpoint of interest is liver damage, as measured by serum sorbitol dehydrogenase (SDH) 
levels in female CD-I mice after 14 days of exposure by oral gavage. Seven separate 
experiments were conducted: four single chemical studies and three mixture studies. A 
more complete description of the experimental methods is provided in Gennings, et al. 
(1997; 2000). The mixture studies were all conducted with a fixed ratio of the four 
10 chemicals (a ratio of (0.34: 0.29:0.32:0.05) for (BDCM:CDBM:CHCl 3 :CHBr 3 )) with 
three different total dose values (0.05, 1 .5, 3.0 mM/kg) and appropriate positive controls 
(a fixed dose of each chemical alone for quality control purposes). The fixed mixing ratio 
of interest is associated with the chlorination process of disinfection (McDonald et al., 

1 999) . The purpose of this study was to determine if the four DBPs combine additively, 
1 5 i.e., with no interaction. 

In Section 2, we provide further details about the DBP study and results. Section 
3 develops the model along a fixed ratio ray in the generalized linear model framework 
using quasi-likelihood estimation. In addition, tests for the null hypothesis of additivity 
are developed and described. Experimental designs involving mixtures of 

20 drugs/chemicals along fixed ratio rays using D-optimality and D s -optimality criteria are 
presented in Section 4. Two-stage designs are also described in Section 4. We apply the 
methodology developed throughout Sections 3 and 4 to the DBP data and present the 
results in Section 5. Finally, in Section 6 we conclude with a brief discussion. 
2. Example Description 

25 The data used to illustrate the methodology developed in the following sections involves 
the combination of four chemicals, which are disinfection by-products in drinking water. 
Summary statistics for the single chemical data are provided in Gennings et al. (1997, 

2000) . A generalized linear model with a power link function, 

g(/Am, add)= ML = Po + E yi X i ' (1 ) 
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where c=4 is the number of chemicals, jc, is the dose of the / th chemical, (3q is the 
unknown intercept, # is the unknown slope parameter for the /' th chemical, was fit to the 
single chemical data to produce an additivity model for the SDH response. Following 
Gennings et al. (1997; 2000) the variance of SDH was assumed to be related to the mean 
5 (i.e., tm for t a dispersion parameter), and the transformation parameter in the power link 
was set to 0.5 for the analysis of these data. The resulting parameter estimates (using a 
quasi-likelihood estimation criterion) and associated p-values are provided in Table 1 . 
All four chemicals are associated with an increase in SDH as the dose of the chemical 
increases (Figure 21). In addition, the additivity model seems to adequately represent the 
1 0 single chemical data as evidenced in Figure 2 1 . 



Table 1: Predicted model parameters for the additivity model given in (1) using the 
single chemical data only. The transformation parameter, X, was set to 0.5 as used in 
Gennings etal. (1997). 



Parameter 


Estimate 


Standard Error 


P-value 


Po 


4.402 


0.213 


O.001 


Yi (BDCM) 


3.567 


0.282 


O.001 


y 2 (CDBM) 


3.833 


0.316 


O.001 


Y3 (CHCI3) 


1.538 


0.239 


O.001 


Y 4 (CHBr 3 ) 


2.441 


0.231 


O.001 


t (dispersion parameter) 


5.293 







15 

After the single chemical studies were completed, three mixture studies were 
conducted. These consisted of mixture points, vehicle controls, and positive controls 
(each chemical alone at 1.52 mM/kg for BDCM, CHC1 3 , and CHBr 3 ; 0.76 mM/kg for 
CDBM). In order to verify that the dose-response relationship was consistent across the 

20 different studies, the additivity model given in (1) was used to predict the response for 
each of the positive control groups. Table 2 provides comparisons of the observed 
sample means at each of the positive control groups and the predicted response using the 
additivity model in (1). The p-value is associated with the Wald-type test for equivalence 
of the two means. Using Hochberg's correction for multiple comparisons, six of the 

25 twelve comparisons resulted in the conclusion that there was a significant shift in the 
dose-response relationship between the original single chemical dose-response studies 
and the mixture studies. Therefore, it was concluded that the single chemical data should 
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not be used to describe the additive or 'zero interaction' case along the mixture ray of 
interest. 



Table 2: Analysis of Positive Control values: Observed sample means from the positive 
5 controls in the mixture study and predicted mean responses from the single chemical data 
in the additivity model as defined in Table 1 . 



Chemical 


Dose (mM/kg) 


Mean SDH 
(STD) 


Predicted from 
Single 

Chemical Data 
(SE) 


P-value 


mixture 
study=0.05 
total dose 










BDCM (n=9) 


1.52 


107.3 (97.5) 


96.5 (7.80) 


0.747 


CDBM (n=9) 


0.76 


38.0(11.2) 


53.5 (3.57) 


0.004* 


CHC1 3 (n=8) 


1.52 


63.1 (34.9) 


45.4 (4.46) 


0.182 


CHBr 3 (n=7) 


1.52 


64.0 (27.8) 


65.8 (5.20) 


0.878 


mixture 
study=l .5 total 
dose 










BDCM (n=9) 


1.52 


56.0 (17.4) 


96.5 (7.80) 


O.001* 


CDBM (n=8) 


0.76 


28.1 (6.7) 


53.5 (3.57) 


O.001* 


CHCI3 (n=6) 


1.52 


46.6 (21.5) 


45.4 (4.46) 


0.905 


CHBr 3 (n=8) 


1.52 


56.0 (24.8) 


65.8 (5.20) 


0.340 


mixture 
study=3.0 total 
dose 










BDCM (n=6) 


1.52 


48.8 (23.3) 


96.5 (7.80) 


0.0002* 


CDBM (n=7) 


0.76 


19.3 (3.6) 


53.5 (3.57) 


O.001* 


CHCI3 (n=8) 


1.52 


24.3 (5.8) 


45.4 (4.46) 


O.001* 


CHBr 3 (n=8) 


1.52 


86.4 (89.2) 


65.8 (5.20) 


0.521 



* Significant at the 5% level using Hocherg's correction for multiple testing. 



The objective of this paper is to develop a method for testing for departure from 
10 additivity when the single chemical data cannot be used to describe the no interaction 
case. The mixture data of interest are assumed to result from a fixed-ratio ray design. 
Table 3 provides summary statistics for the DBP mixture data described above. The 
three 'total dose' groups represent the mixture of the four chemicals at a fixed ratio of 
(0.34: 0.29: 0.32: 0.05) for (BDCM: CDBM: CHC13: CHBr 3 ), but with total dose values 
15 of 0.5, 1.5, and 3.0 mM/kg. 
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Table 3; Summary statistics for the mixture data. 



Mixture group 


Total dose 


Mean SDH 


STD SDH 


Sample size 


Control 


0 


17.42 


3.43 


27 


Chlorination* 


0.05 


19.61 


4.36 


17 


Chlorination* 


1.5 


46.29 


19.26 


17 


Chlorination* 


3.0 


76.22 


12.10 


16 



* Mixing ratio of (BDCM: CDBM: CHC1 3 : CHBr 3 ) was (0.34:0.29:0.32:0.05). 



3. Model and Estimation 

5 Suppose we are interested in studying the interaction among c drugs/chemicals in 

combination at a fixed ratio. Let y ijk be the kth observation of the jth dose of the zth 
chemical; 1 < i < c, 1 <j < d t , 1 < k< n,y; and N=^. . riy . Denote the mean of Y as 

E(Y)= ju and the variance of Y as Var(Y)=TV(u), where V(u) is an assumed known 
function of the mean and x is an unknown scale parameter. Using a generalized linear 
10 model, the mean of Y is associated with the doses of the c chemicals through a link 

function (McCullagh and Nelder, 1989). With the assumption of the form for the mean 
and variance of Y, quasi-likelihood methods can be used to estimate the model (mean) 

parameters (e.g., McCullagh, 1983), i.e., Q(u,r;y)= f y ~ W d w is maximized across 

* tV{w) 

the N observations in the study using an iterative algorithm. A method of moments one- 
15 step estimator is used for t (e.g.,. McCullagh and Nelder, 1989). 
3.1 Model along a Fixed Ratio Ray 

For a combination of c drugs/chemicals with mixing ratio [ay: a?.. . .:a c ], the additivity 
model can be expressed as given in (1) with the additional assumption that var(Y) = 
tV(u). If the model is parameterized to include the terms associated with all 2, 3,. . .,c- 
20 drug/chemical interactions, then the interaction model assumes the following form 

c c c c c c 

^^interaction) = A) + £ M + S ZIW* 

i = l i = \j = \ i = \j = \k = \ 

i<j i<j<k 

c c c c 

+ S Z Z E Yijkl*iXjXkXl+ - + Y\2...c x \ x i- X c ( 2 ) 
,=1 j=\k=\l=\ 

i<j<k<l 

where x, is the dose of the ith drug/chemical, (5$ is the unknown intercept and # is the 
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unknown slope parameter for the z'th drug/chemical, yij is the unknown parameter 
associated with the interaction between drugs/chemicals i and j, is the unknown 
parameter associated with the interaction between drugs/chemicals i,j and k, and yn...c is 
the unknown parameter associated with the c-factor interaction among the 
5 drugs/chemicals. Pure higher order terms and associated cross-product terms are not 
included in (2) as such terms do not allow the corresponding additivity model (i.e., the 
model in (2) where all cross product terms are zero) to satisfy Berenbaum's definition of 
additivity. 

When the combination data are collected along a fixed ratio ray with mixing ratio 
10 {ai:ar. ... :a c ), then a mixture point [xi, X2, ,x c ] is uniquely defined by [ait, a,2t, ... a c t] 

where t is the total dose. The model of no interaction defined in equation (1) reduces to 
the linear predictor being in the form of a simple linear model in terms of total dose, i.e., 

&(Pm,add) = A) + Pit, (3) 

where j}\ = ^Ty ,-a,. . Additionally, the model that allows for interaction, defined in 

1 5 equation (2), can be expressed with the linear predictor portion of the model being a 
higher order polynomial in terms of the total dose, defined by 

^^interaction) = A) + A* + Z A*'. < 4 > 

;=2 

where A = Z fr a / > & = Z Z TtyV; > A = Z Z Z rijk a i a j a k > etc. Here, fa is 

/=1 <=1 j=\ i=l j=lk=\ 

i<j i<j<k 

the unknown parameter associated with all two-factor interactions, fa is the unknown 
20 parameter associated with all three- factor interactions, . . ., and fa is the unknown 

parameter associated with the c-factor interaction. With mixture data only along a fixed- 
ratio ray, the entire response surface cannot be adequately (if at all) estimated. Thus, the 
g parameters in (2) are not estimated. Instead, the design supports the estimation of the 
polynomial parameters in (4). 
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3.1 Inference 

Let 0 =[fio, P\, Pi>. . ., flcY, where 0 is found by replacing unknown parameters 
with quasi-likelihood parameter estimates. Following McCullagh (1983), the large 
5 sample covariance matrix of 6 is var(#) = [I(0)]~ l , where the expected quasi- 
information matrix evaluated at 9\s written as 

I(^) = D'(rV(//))- 1 D|. (5) 

where 









dPo 
























dfio 







"Vfa) o ... 0 

0 V(^ 2 ) ••• 0 
10 V(//) = : ; : • 

0 0 - V(^) 

It is important to note that the information matrix in (5) depends on the design of the 
study (in terms of dose location and subject allocation) through the rows of D. This is 
exploited in section 4 to determine optimal experimental designs. 

Of ultimate interest is the determination and classification of departures from 

15 additivity among a particular combination/mixture of chemicals. When comparing the 
additivity model to the interaction model, defined in (3) and (4) respectively, evidence of 
curvature on the link scale indicates departure from additivity; i.e., there is interaction 
among the chemicals if at least one #*0, i=2,...,c. When 0=[6\ : 6z]', where 6\ =[j%, 
P\]' and 02=[P2,-- ■, Pc]', the test for departure from additivity can be written as 

20 Ho: ft = 0vs.Hi: ft*0. (6) 

In testing such an hypothesis, the analyst generally selects from three large sample tests: 
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the quasi-likelihood ratio test (McCullagh, 1983), the score test, and the Wald test (e.g., 
McCullagh andNelder, 1989). 

4. Experimental Designs 

5 We have shown that when a ray design is used the hypothesis of additivity is 

equivalent to the hypothesis of the adequacy of the model where the linear predictor is in 
the form of a simple linear model. In other words, the test for additivity is testing that the 
parameters associated with the higher order polynomial terms in the linear predictor are 
equal to zero. The tests of the hypothesis in (6) either implicitly or explicitly involve the 
10 significance of 0 2 terms relative to Var(# 2 ) • Thus, designs which lead to precise 
estimation of 0 2 preferred. Hence, a reasonable criterion for use when designing 
mixture studies to detect interactions requires optimizing the precision of the estimates of 
the model parameters. Such designs are associated with an increase in the power of the 
test of additivity. 

15 

4.1 D-Optimal and Ds-Optimal Designs 

Define t=[t 1 t 2 ... t c+r { as the (c+r)x 1 vector of total dose design points 
where the number of total dose levels is fixed. The minimum number of levels needed to 
test the hypothesis of additivity is (c+1). Define q = q 2 ... q c + r \ as the vector of 

20 sample size allocations where ]T# f =1 and the number of observations in the i th dose 

group is nj=N<jrj. The D-optimal design is the design that minimizes with respect to t and q 
the generalized variance of the estimated model coefficients (Kiefer and Wolfowitz, 
1959), where the generalized variance is defined as the determinant of the variance- 
covariance matrix. In the present context, the D-optimal design would maximize the 

25 determinant of the Fisher's quasi-information matrix; i.e., 

, (?) 

T I 

where I(/7) is the Fisher's quasi-information matrix defined in equation (5). 
Maximization of the determinant of the information matrix is accomplished by adjusting 



max 
t,q| 
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the rows in D in (5) for the total dose values and the number of observations associated 
with each dose group at each step in a direct search algorithm. 

In Section 3 we introduced methodology for testing the hypothesis of additivity 
given a generalized linear model framework, using quasi-likelihood estimation. We 
5 showed that the hypothesis of additivity is rejected when higher order polynomial terms 
are required in the total dose-response model. Thus, it is important that we have precise 
parameter estimates for these higher order polynomial terms in the linear predictor. We 
will develop methodology for finding a D s -optimal design based on this subset of terms 
in the linear predictor. 

1 0 The Ds-optimal design is used to select design points that minimize with respect 

to t the generalized variance of the parameters associated with the higher order 
polynomial parameter terms in the linear predictor; that is, 

minVar(0 2 ), (8) 

where Var(0 2 ) is the appropriate subset of the inverse of I(0\, fy), defined in (5), and is 
15 evaluated at specified values of &\ and (h- In practice, situations may arise where the 
search needs to be constrained to a particular region of interest. In either case, the 
Nelder-Mead direct search algorithm (Nelder and Mead, 1965) can be used to accomplish 
the minimization procedure. 

The approach discussed above results in a D s -optimal design. However, the 
20 practical use of a Ds-optimal design may be questioned since the design is dependent 

upon unknown parameters. Thus, the scientist who employs the optimal design may need 
to "guess" the parameter values. The "guess" may be based on preliminary studies or 
information found in the literature. Another approach is to develop two-stage designs 
that provide good statistical properties, as suggested by several researchers, including 
25 Abdelbasit and Plackett (1983), Minkin (1987) and Myers et al. (1996). 



4.2 Two-stage Designs 

A two-stage design may be employed when there is only scant information on 
which to base a reasonable "guess" of parameter values and when the researcher believes 
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the experimental conditions in two studies are reproducible. A two-stage procedure uses 
the second stage to complement the first stage. 

Following the approach taken by Minkin (1987), the total quasi-information 
matrix of the entire experiment can be expressed as the sum of the information matrices 
5 of the first stage and the second stage given the first stage, i.e., 



Minkin (1987) based this result on conditional likelihood theory; i.e., the fact that the log 
10 likelihood of a two-stage design can be expressed as the sum of the log likelihood of the 
first stage and the log likelihood of the second stage given the first stage. 
The following steps outline the two-stage design. 



Vai = A, P c ) + — I.I/CA. A. A.-. A). 




(9) 
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(1) Provide initial parameter estimates based on previous experimentation, 
a pilot experiment conducted specifically for this purpose, or a guess 
by the scientist (J3 0{ , fa, fhf,..., Pcf )■ 



(2) Choose the total sample size for the entire experiment, N= Nf+ N s . (Nf 
is the number of experimental units allotted for the first stage and Ns is 
the number of experimental units allotted for the second stage). 



(3) Based on Nf, Pot, flu, (hu---, Pa , apply, for example, the D-optimality 
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criterion defined in equation (9) to obtain the values of total dose and 
sample size allocation (i.e., q) that maximize the determinant of the 
quasi-information matrix. This is accomplished by using the Nelder- 
Mead direct search algorithm. 



25 



(4) Perform the first stage of the experiment at the first stage D-optimal 
total dose values (?i f , . .,^ r i)f) with Ni(# If , q 2 u. . .,3Vri)f) 
observations at each dose group and obtain maximum quasi -likelihood 
estimates for the model parameters (/? of , A f , J3 2 f ,—, P c {) . Equal 



allocation in the first stage may be considered reasonable. 
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(5) 



Based on the first stage parameter estimates (J3 0{ , /? lf , /9 2 f P c f )> 



if 9 ^s, Ns, apply the Ds-optimality criterion defined in equation (8), 



5 



where Var(0 2 ) is replaced with the appropriate subset of the inverse 
of the total conditional information for the two-stages, I t0 t 3 i, given in 
equation (9), to obtain the second stage total dose values (*i S |f, tj^,. . ., 

*(c+r2)s|f). 



(6) 



Conduct the second stage experiment at the second stage Ds-optimal 
total dose values (t\ s \f, fc S |f,. • •, J(c+r2)s|f) with ni=N2<?i observations at the 
i th dose group, i=l,.. .,c+r2. 



10 



(7) 



Based on the observed responses, y iJk , at tu, ^ljf, ti& f 2s|f,..., 
^r2)s|f, estimate the final model parameters 0 X and 0 2 using maximum 



quasi-likelihood estimation. In addition, the estimate of Ts|f is found 
by using moment estimation. 
This will result in a two-stage D-Ds-optimal design consisting of 2*c+rl+r2 
15 design points. 

5. Example 

Descriptive statistics for the mixture of the four DBPs described in Section 2 are 
provided in Table 3. Clearly, as the total dose increases there appears to be an increase in 

20 the mean SDH. In addition, the standard deviation appears larger for the total dose 
groups associated with a higher mean response. Thus a quasi-likelihood criterion for 
parameter estimation with Var(Y)=7yi seems reasonable for these data. Using the same 
link function as considered in the additivity model (i.e., g{\i)= /i X where X=0.5), a 
quadratic model was initially fit to the data (Table 4a, Figure 22). From (3), the 

25 additivity model along the fixed-ratio ray includes an intercept and slope parameter. 

Assuming the underlying dose-response surface associated with the four chemicals is as 
described in (4), the test for the significance of the quadratic term along the fixed-ratio 
ray is analogous to the simultaneous test of all pairwise interaction terms. The p-value 
associated with the significance of f} 2 (Table 4a) is 0. 1 52. Although not significant at the 

30 usual 5% significance level, the p-value is compelling enough to motivate further 
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investigations of the mixture. With a larger design, it would be possible to again test for 
the significance of all higher order terms, which is associated with the test of additivity. 

Table 4a: Estimated model parameters using the mixture model given in (2) for a 
5 quadratic generalized linear model with power link with the transformation parameter 



X=0.5. 



Parameter 


Estimate 


Standard Error 


P-value 


Po 


4.234 


0.124 


O.001 


P. 


1.938 


0.300 


O.001 




-0.146 


0.102 


0.152 


T 


2.58 







Table 4b: Estimated model parameters using the mixture model given in (2) for a cubic 



Parameter 


Estimate 


Standard Error 


P-value 


Po 


4.174 


0.155 


O.001 


Pi 


5.246 


5.250 


0.318 


P2 


-3.415 


5.179 


0.510 


P3 


0.724 


1.147 


0.528 


T 


2.60 







* The likelihood ratio simultaneous test for the significance of P2 and P3 was not rejected 
(p=0.294). 



The design used for the mixture data (Table 3) was selected to include a small 
15 dose value, a large dose value and somewhat of a midpoint value. It was not based on a 
formal design strategy. Therefore, it is of interest to propose a design that might be 
useful in further investigations of the mixture. Two options are available for 
consideration. First, the observed data to date could be considered preliminary data and a 
full new study could be conducted. The parameter estimates found in the observed data 
20 could be used in a one-stage D s -optimal design. In this approach, the new data would be 
analyzed separately from the original mixture data. Another approach is to consider the 
observed data as a first-stage and patch these data with additional dose groups so that a 
fourth degree model could be fit to the combined data. A fourth degree model is of 
interest as the mixture includes four DBPs and such a model allows for a test of up to a 
25 four-way interaction among the chemicals along the fixed-ratio ray of interest. Of 
concern is that the dose-response relationship shifted from the initial single chemical 
studies to the mixture studies. If another shift occurs between the first and second stages, 
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then the data could not be combined for a full analysis. Both approaches are considered 
in more detail below. 
One-stage design 

A Ds-optimal design associated with the test for interaction is a design that 
minimizes the volume of the confidence ellipsoid about the higher order parameters. In 
constructing this design, the mixture data were refit with a cubic model with the same 
power link function (Table 4b). Assuming an underlying dose-response surface as given 
in (4), the test for interaction would be based on the significance of ft and ft (i.e., the 
quadratic and cubic terms which are associated with all pairwise and three-way 
interactions). The design algorithm resulted in values for the total dose groups and 
sample size allocations associated with minimizing the determinant of the 2x2 covariance 
matrix associated with the higher order terms. A region constraint was added to the 
Nelder-Mead algorithm that required the largest total-dose group to be no larger than 3.0 
mM/kg. This was selected based on the range of the single chemical data and because 
the investigator is more interested in the low-dose region than a high-dose region. The 
resulting design is provided in Table 5. 



20 



Table 5: One and two-stage designs associated with testing for interaction among the 
four DBPs using the model parameters in Table 4b to describe the assumed shape of the 



Design 


k d x d 2 d 3 dA 

[<lo 9l <l2 ?3 94 J 


One stage design 
(assuming vehicle + 
four points) 


f 0 0.9 1.8 2.5 3.0l „ VT , 

^ > for N observations 

(0.26 0.14 0.16 0.21 0.23J 


Two-stage design 
(assuming observed 
data is first stage and 
adding new vehicle+2 
points; N=Ni+N 2 ) 


f 0 0.05 1.5 3.0l ^ , 
First stage: < > for Ni observations 
[0.35 0.22 0.22 0.21J 

f 0 1.0 1.9] 
Second stage: < > for N2 observations 
5 [0.57 0.28 0.15J 
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Two-stage design 

A two-stage design assumed the observed mixture data was the first stage. The 
idea is to 'patch' this design with two additional points, so that a quartic model could be 
estimated. With four chemicals in the mixture it is of interest to test for the significance 
5 of up to a four-way interaction. Assuming the underlying dose-response surface as 
defined in (2), a four-way interaction is associated with a fourth degree term along the 
fixed-ratio ray. 

The two-stage design is based on the conditional information matrix given in (9). 
The first stage information is added to the conditional second stage information matrix. 

10 This was accomplished assuming the underlying parameters are as given in Table 4b and 
the first stage information was the inverse of the covariance matrix of the model 
parameters. Again, a region constraint was included in the Nelder Mead algorithm 
requiring the largest total dose was no greater than 3.0mM.kg. The resulting second- 
stage design is provided in Table 5. In a two-stage design, both the first and second-stage 

1 5 data are combined in the final analysis. Of course, preliminary checks for the equality of 
corresponding positive control group means would be conducted. 
6. Discussion 

A goal of this paper is the development of methodology for detecting departures 
from additivity among mixtures of drugs/chemicals taken along a fixed mixing ratio 

20 when a generalized linear model using quasi-likelihood estimation is used. This 
approach is useful when interest is focused on a fixed ratio ray, when resources are 
limited, and when only the first two moments of the distribution of the response variable 
can be specified. Importantly, the methodology developed in this paper does not require 
single-drug/chemical data to test for departures from additivity along the fixed ratio ray 

25 of the drugs/chemicals. The present example illustrates a situation where single chemical 
data were available but were unusable due to a shift in the dose-response curves. While 
single chemical data would be desired they are not necessary. 

The number of datapoints necessary to test for interaction along a fixed-ratio ray depends 
on the degree of the polynomial assumed. With c chemicals in the mixture, one may be 
30 interested in fitting a c-degree polynomial. The minimum number of dose points would 
be c+\. However, as the number of chemicals in the mixture increases, the analyst may 
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be willing to assume that any higher-order interaction greater than say of degree k is 
negligible. In this case the minimum number of data points is k+1. Of course, additional 
data points can be added to further support the shape of the dose-response relationship. 
An alternative approach to test for chemical interactions along a fixed ratio when 
5 single drug/chemical data are available is described in Gennings et al. (2002). These 

authors compared the estimated slope of an additivity model along a fixed mixing ratio of 
interest to the estimated slope based on combination data along the ray. Here, the single 
chemical dose-response data are experimentally observed to support the estimation of an 
additivity surface. However, as the number of drugs/chemicals increases, the likelihood 
1 0 that the single drug/chemical data are available decreases. 
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EXAMPLE 8. Analysis of an Interaction Threshold in a Mixture of Drugs and/or 
Chemicals 



1. Introduction 

5 A major task when studying chemical mixtures is to determine whether departures 

from additivity, i.e. interactions, among the chemicals in the mixture exist. Current 
methodology (Gessner and Cabana, 1970, Carter, et.al, 1988, Kelly and Rice, 1990, 
Gennings, 1995, Dawson, Carter, and Gennings, 2000) provides foundations to study and 
characterize interactions by utilizing concepts involving isobolograms, statistical models, 

10 and the interaction index (Berenbaum, 1981). Toxicological data may suggest dose- 
dependent interactions (e.g. Konemann and Pieters, 1996, Gennings, et.al. 2002) where 
the existence and nature of the interaction may change with dose. The U.S. EPA (1996, 
2000, 2002) and Carpy, Kobel, and Doe (2000) suggested that low-dose regions of 
mixtures of chemicals should be associated with additivity, while interactions might 

1 5 occur in higher-dose regions. 

The definition of additivity used here is consistent with the classical isobologram 
(Loewe and Muischnek, 1926; Loewe, 1953). For a mixture of two agents, an 
isobologram is a contour of constant response of the dose-response surface superimposed 
on the line of additivity defined by the equi-effective levels of the individual components 

20 of the mixture. When the observed contour is below the line of additivity, a synergism 
can be claimed. When the observed contour is above the line of additivity, an 
antagonism can be claimed. When the observed contour is coincident with the line of 
additivity, additivity can be claimed. For a mixture of c > 3 chemicals, the production of 
an isobologram suffers from the difficulties associated with displaying c-dimensional 

25 figures. 

For a mixture of c chemicals, an additivity model can be written as 
g{n) = /3 0 + fi x x x + J3 2 x 2 +... + /3 c x c 
where xj, x 2 , .... x c are the doses of the c individual chemicals, g(^i) is a user-specified 
link function (McCullagh and Nelder, 1989), and /? 0 , y# c are unknown parameters. 
30 At a fixed response, uo, g(iu 0 )-fi 0 = P\* x + P 2 X + • • • + P c m & 
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1 = 



X, 



x c ^_ x x 




= Interaction Index 



ft 



+ ■■• + 



g(Mo)-Po ED® 

ft 



(1.1) 



This follows since 



g(Vo)-ft 
ft 



= ED™ , the dose associated with the response uo for the 



i th component of the mixture. Berenbaum (1981) related the interaction index to the 
5 isobologram and showed that when the interaction index is not equal to one, an 
interaction, i.e. departure from additivity, is present. The interaction index is important 
since it does not have the graphical limitations associated with producing plots of multi- 
dimensional isobolograms when the number of chemicals in the mixture is .greater than 
three. 

10 Carter, et.al, (1988) showed that when model parameters associated with 

interaction in a generalized linear model are different from zero, the interaction index is 
different from one. For example, consider a two-chemical mixture with interaction, i.e. 



When Pi 2=0, the interaction index is equal to one. For g(/j)> ft 0 , if the estimate of 
Pi2>0, the interaction index is less than one which is indicative of a synergism. Similarly, 
when pi2<0, an antagonistic interaction can be claimed. 

Gessner and Cabana (1970) studied the mixture of ethanol and chloral hydrate, 

20 using the loss of righting reflex in mice as the response. These authors compared 
confidence interval estimates about experimentally determined EDso's at nineteen 
different dose combinations to lines formed by connecting the lower confidence limits of 
the ED 50 's of the individual drugs. They interpreted these latter lines as some form of a 
confidence bound about the line of additivity formed by connecting the ED 50 's 

25 determined for each drug. Evidence for a departure from additivity was provided 
whenever the confidence bounds for the experimentally determined ED 5 o for a 



g(M) = Po + M + A*2 + #2*1*2 



From this, it follows that, 



15 



1- 
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combination did not overlap with the "confidence" region for the line of additivity, as 
shown in Figure 1 (taken from Figure 1 of Gessner and Cabana, 1970). 

Although the statistical approach used to analyze their data is crude, their work is 
important because it provides an empirical estimate of the ED 5 o contour of the underlying 
5 dose response surface associated with these two chemicals. It is noteworthy that the 
experiment used a total of 1681 animals in 21 treatment groups, the 19 combinations of 
the two chemicals mentioned previously, and two single chemical groups (Gessner and 
Cabana, 1967). 

The isobologram of Gessner and Cabana (Figure 1) associated with the empirical 
10 EDso's of ethanol and chloral hydrate is consistent with the existence of an interaction 
threshold. Note that the mixture appears to be additive at the 50% level of response for 
chloral hydrate less than 125 mg/kg and ethanol greater than 1200 mg/kg. When chloral 
hydrate exceeds 125 mg/kg and ethanol is less than 1200 mg/kg, there is a synergistic 
interaction. This suggests the possibility of an interaction threshold when the drugs are 
15 combined. 

The boundary that separates the dose space into regions of interaction and 
additivity is of interest to locate. The development and subsequent analysis of a model 
that accommodates the elucidation of the boundary of the two regions is the objective of 
this report. The interaction threshold boundary may take a variety of different shapes, and 

20 the shape of this boundary is not likely to be known. Our goal is to develop a general 
procedure to accommodate various potential shapes for this boundary. 

In Section 2, we develop a dose-response model that allows an interaction 
threshold boundary. In Section 3, we describe parameter estimation procedure and tests 
of hypotheses of specific interest. An example illustrating these methods using a mixture 

25 of ethanol and chloral hydrate is presented in Section 4. A discussion of these methods 
and some concluding remarks are presented in Section 5. 
2. Development of the Interaction Threshold Model on the Response Surface 

Suppose that we are interested in estimating an interaction threshold boundary for 
a mixture of c chemicals where both single chemical and mixture data are observed. Let 

30 Vjjk be the k th observation of the j th dose of the i th chemical and let the mean of Y be 
E(Y)=p. Assume that var(Y)=xV((j,), where V(u) is assumed to be a known function of 
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the mean and x is an unknown scale parameter. Given a c-dimensional combination, we 
assume that the interaction threshold boundary can be defined such that the value of the 
c th component can be expressed as a function of the remaining c-l chemicals, i.e. 

x c =Q(x„x 2 ,...,x c _ l ). 
2.1 The General Interaction Threshold Model 

Consider the following generalized linear interaction threshold model describing 
the relationship between the response and the doses of the c chemicals in combination: 



r=l I 
c c-l c c-2 c-l c | 

A +2X*r +Z5X*,*, + ZZ ftw, +• • -Pi.^h ..x c x c >g < \x 2 ,. . ^ 

k r=l r=l s=r+l r=l s=r+l<=3+l J 

(2.1.1) where xi, X2, x c are the doses of the individual chemicals, g(u) is a user- 
specified link function (McCullagh and Nelder, 1989), and /?„,#,. ..ytf c ,# 2 ,...,y9 12 c are 
unknown parameters. The model can be made continuous by requiring the values of g(u) 
at the threshold boundary to be equal, i.e., at the boundary 

\r=\ VVs=r+l J J r=ls=r+l \\t=s+\ J J J 

(2.1.2) 

For simplicity, consider a mixture of two drugs/chemicals. In this case, the 
continuity constraint given in (2.1.2) requires that fi n x x Q{x x ) = 0 . Thus, for the two- 
drug/chemical case, (2.1.1) can be written as 

goo = & + /?,*, + p 2 x 2 +|Ai(*,(*2 -0(*i)))K(* 2 

(2.1.3) 
where -£>(*,)) ^ 

2.2 Several Formulations of a Threshold Region When c=2 

It is unlikely that investigators will know the form of Q(x\) or, equivalently, the 
location of interaction thresholds, if they exist. When the form of Q(x\) is specified in the 
interaction threshold model, it is possible to estimate the unknown constants that 
parameterize the interaction threshold boundary. Four different functional relationships 
are considered here and will be used in the example given in Section 4. 
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(i) Segmented Line (Figure 2 A) 

For this boundary, Q\(xi) is given by 



^ / % f A + Bx x if x t < C 

^--{d + E;, if J|2 C (221) 
with the continuity constraint that A + BC = D + EC . After application of the method of 
Gallant and Fuller (1973), (2.2.1) becomes 

Q x (x x ) = A + Bx x + [(E - B)(x, - C)]/ + (*, - C) (2.2.2) 



where 7 + U-C) = ^ 



1 0 (ii) Nonlinear boundaries 

(a) Four parameter logistic boundary (Figure 2B) 



Here, 

15 Q 2 {x x ) = A + Bf{x x ) (2.2.3) 

For this example, we have chosen fix, ) = , the logistic CDF. Any 

cumulative distribution function would suffice to impose the boundary we envision. This 
boundary imposes asymptotic boundaries at A and A+B, for the levels of the second 
drug/chemical associated with the threshold boundary. 
20 (b) Inverse Cubic Boundary (Figure 2C) 

Here, 

QJx.) = A + l - — (2.2.4) 

31 B(x l -E) + C(x i -E) 2 +D(x l -E) 3 

This relationship permits asymptotes to minimum values of each of the two 
25 drugs/chemicals, E for drug/chemical 1 and A for drug/chemical 2. 
(Hi) Elliptical Boundary (Figure 2D) 
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It may be the case that for a combination of two drugs/chemicals, interaction may 
be constrained to a region within the dose space. Consider the situation where this region 
is elliptical. Here, if Ax 2 +Bx 2 2 + Cx { +Dx 2 + Ex x x 2 +1<0, interaction is present. 
When Ax 2 + Bx 2 2 + Cx, + Dx 2 + Ex x x 2 + 1 > 0, there is no interaction. That is, when the 

5 combination of doses, (xi,X2), is inside the ellipse, the combination is associated with 
interaction. When the combination of doses, (xi,X2), is outside the ellipse, the 
combination is associated with no interaction. 

The nature of this region of interaction requires that the user modify the general 
interaction threshold model to preserve continuity. For this situation, we propose the 

10 following continuous model 

gOu) =B 0 + B x x x + B 2 x 2 + [/ViV]'-(04(*i.* 2 )) ( 2 - 2 - 5 ) 
where Q 4 (x, , x 2 ) = Ax 2 + Bx 2 + Gt, + Dx 2 + Ex x x 2 + 1 and 
[0 Q 4 (x„x 2 )>0 



The values of the variables x, and x 2 are scaled such that the relationship is continuous 
15 on the interaction threshold boundary. There are multiple ways to accomplish this, and 
we discuss one such method in Section 4. 

For the interaction threshold model that has been defined for the various 
formulations of Q, our objective now is to estimate the model parameters. 
3. Estimation and Hypothesis Testing 
20 3.1 Parameter Estimation 

For the combination of two chemicals, the interaction threshold models are 

g(M) = A> + A*. + 02*2 + [A 2 (*.(*2 " <2(*,)))]/ + (* 2 - 
(3.1.1) 

for Q(x,) given by (2.2.2), (2.2.3), and (2.2.4), and 

25 g(M) = 0 O + A*, + P 2 x 2 + [Az*. V]'- (GO,,**)) (3-1-2) 

for Q(x,,x 2 ) given by (2.2.5). 

The methods of maximum quasi-likelihood (Wedderburn, 1974; McCullagh and 
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Nelder, 1989) can be used to estimate the unknown parameters from the interaction 
threshold model. The log-quasi-likelihood function is 

\ vV{t) 

Assuming independence across observations, the total log quasi-likelihood is 

q(/i,t;y) = , t q(My >r;y ijk ). (3.1.3) 

A direct search algorithm, such as the Nelder-Mead algorithm (1965) can be used 
to estimate the unknown parameters due to its ease with which it can incorporate the 
constraints associated with the interaction threshold boundary. We use a closed form 
estimator for x suggested by McCullagh and Nelder (1989) shown as 

(y uk -^ = _X^_ 
6iV(M 9 )(N-p) N-p 



where N is the total number of observations, and p is the total number of parameters 
estimated in the particular model of interest. X 2 is the generalized Pearson statistic 
asymptotically distributed % 2 with N-p degrees of freedom (McCullagh and Nelder, 

15 1989). Estimation is performed in SAS, Version 8.2 using Proc NLP specifying 
maximization of the quasi-likelihood by the Nelder-Mead algorithm, and Proc ML for 
estimation of f and the appropriate variance-covariance matrix. 
3.2 Hypothesis Testing 

Separate tests of the hypotheses of overall additivity and the hypothesis of the 

20 presence of the interaction threshold boundary are performed. These hypotheses are 
tested using a log quasi-likelihood ratio test (McCullagh and Nelder, 1989). Denote the 
maximum values of the log quasi-likelihood under H 0 and H A as 
Q(y;0,r ) and Q(y,0, f) , respectively. Then, under Ho, 

QLR = 2^(y;0) - Q(y\6)\ is approximately distributed x%\ (McCullagh, 1983), where 
25 k is the difference in the number of parameters between the restricted and unrestricted 
models. Following McCullagh and Nelder, 1989, f is a consistent estimator of the 
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unknown parameter t and is approximately distributed — — . 

n-p 

Hinkley (1969) recommended the use of the likelihood-ratio statistic for 
inferences about the join point or joint inferences about the join point and model 
parameters in segmented regression. Feder (1975) noted that the distribution of the 
5 likelihood ratio might be better approximated in finite samples by a multiple of an 
asymptotically equivalent F-distribution (Seber and Wild, 1989, pgs. 451-452). 
OLR 

Following this logic, we use — — , where f is calculated in the unrestricted model, as a 
ki 

likelihood-ratio test statistic, approximately distributed Fk >n .p. 

(i) The Test for Additivity 

1 0 Partition the interaction threshold model parameter vector from (3.1.1 or 3. 1 .2) as 

© = [©, :0 2 ] where 0, =[/3 0< /3J , 0 2 =[/? 12 ,^], and ^ is the k-vector of model 
parameters associated with the interaction threshold boundary, Q. With respect to a 
mixture of two chemicals, the test of departure from additivity is 
H 0 :© 2 =0 H l :0 2 *0 (3.2.1) 

1 5 For the log quasi-likelihood ratio test, the unrestricted model is given by (3 . 1 . 1 ) or 

(3.1.2) and the restricted model is 

g(M) = A,+A*. +J3 2 x 2 

If the null hypothesis of additivity is rejected, there is evidence of an interaction 
somewhere in the experimental region. It is of interest to determine if an interaction 
20 threshold exists and lies in the experimental region. Testing this hypothesis is equivalent 
to testing if any of the unknown parameters in the interaction threshold boundary, Q, are 
different from zero. 

(ii) The Test for the Presence of an Interaction Threshold Boundary 

We assume that the estimated boundary is within the experimental region. 
25 Partition the interaction threshold model parameter vector as y/ = [^P, : <f> J where 
*P, = [/?„_/?, , P X2 ]' • Tne hypothesis of no interaction threshold is 

#o -1 = 0. H A :l*0 (3.2.2) 
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Again, the unrestricted model is given by (3.1.1) or (3.1.2), and the restricted model is 



If we fail to reject the null hypothesis, it can be claimed that the interaction exists in the 
5 entire dose region. If we reject the null hypothesis, the interaction threshold boundary 
splits the dose region into areas of additivity and interaction. 
4. Ethanol-Chloral Hydrate Example 

Carter, et.al. (1988) provides data from the study of the combination of ethanol 
and chloral hydrate. In their study, the endpoint used was the same as that used by 
10 Gessner and Cabana (1970), the loss of righting reflex. These data, given in Table 1, were 
used to estimate the parameters of the interaction threshold model for each of the four 
considered interaction threshold boundaries. Estimation of parameters in the interaction 
threshold models using the original doses of ethanol and chloral hydrate resulted in a 
wide range of parameter estimates, in some examples, as wide as on the order of 10" 8 to 
15 10 3 . To avoid the appearance of an information matrix ill-conditioned to handle the wide 
range of variances and covariances associated with these parameter estimates, the 
concentrations of ethanol and chloral hydrate are scaled by the conversions 



g(ji) = /? 0 + + 0 2 X 2 + jfl I2 *l*2 




mg 
kg 



chloral hydrate 



425 



which scales the doses to values between 0 



and 1. 



20 



25 



30 
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Since the response variable is binary, g(//) = In — — , and the interaction 
threshold models become 

ln (lf^) = A + A*, + A*2 + LM*i(*2 -0(*,)))fe > 0(*,)) (4-1) 

for gfa) given by (2.2.2), (2.2.3), and (2.2.4), and 

5 = Po+ Px*i + Pi*z + \j3 n x' x x 2 ~\l(Q{x l ,x 2 ) < 0) (4.2) 

for Q(x,,x 2 ) given by (2.2.5). 
To preserve continuity in the interaction threshold model given in (4.2), we 
rescale the doses as 

x x =mm(x l -h 2 (x 2 ),h l (x 2 )-x l ) (4.3) 

10 x 2 =mm(x 2 -v 2 (x x ),v l (x l )- x 2 ) (4.4) 

where 



_ - (Ex 2 + Q + ^j(Ex 2 + Cf - AA{Bx 2 2 + Dx 2 + 1) 
x, - h x (x 2 ) - — 



, - (Ex 2 + C) - J(Ex~ 2 + C) 2 - 4A(Bx 2 2 + Dx 2 + 1) 
x, = h 2 (x 2 ) = — 



/ x - (Ex, +D) + J(Ex x + Df - AB(Ax 2 + Cx x + 1) 
x? = v,(x,) = 

2 IV 1/ ^ 



-(£x. +£>)- J(£x, + £>) 2 - 45(^c, 2 + Cx, + 1) 
15 x, =v,(x,) = - 

2 2V M ^ 



By rescaling the doses in this way, at the interaction threshold boundary 
(boundary of the ellipse), the point is rescaled to zero and thus its contribution to the 



173 



interaction is zero and the model is continuous. Within the ellipse, and x 2 take on 
values that are dependent on the distance from the closest boundary of the ellipse. As the 
point moves closer to the middle of the ellipse from the closest boundary, its contribution 
to the interaction increases, and as it moves away from the middle of the ellipse to the 
5 closest boundary, its contribution to the interaction decreases. 

Each point (xi,X2) is considered individually. The process of rescaling doses 
begins by determining whether the given point lies inside the ellipse. If 
Ax* + Bx 2 2 + Cx x + Dx 2 + Ex x x 2 + 1 > 0 , the point (xi,x 2 ) is outside the ellipse and is 
associated with additivity and is not rescaled. If Ax x 2 + Bx 2 + Cc, + Dx 2 + Ex x x 2 + 1 < 0 , 
10 the point (xi,x 2 ) is inside the ellipse and is rescaled so that its distance from the boundary 
determines its contribution to interaction. An additional constraint imposed on the model 
is that the rescaling of either xi and/or x 2 cannot exceed the value of xi and/or x 2 , 
respectively. This situation occurs if the closest interaction threshold boundary falls in 
negative dose space. 

15 With several different boundary candidates under consideration, the results of 

Pearson Chi-square goodness-of-fit tests are helpful in making choices among the 
candidate threshold boundary forms. Since the p-value associated with the test of fit is 
adjusted for the degrees of freedom, ranking the various forms of Q through the p-values 
may be informative. The results of the goodness of fit tests for each model are contained 

20 in Table 2. 

Table 2 Goodness-of-Fit results for the four considered interaction threshold models 



Interaction Threshold Boundary 
Considered in Model 


Goodness-of-Fit Test 
Statistic 


df 


p-value 


Segmented Line Boundary 


42.33 


31 


0.084 


Nonlinear Logistic-Type Boundary 


36.72 


31 


0.221 


Inverse Cubic Boundary 


44.81 


30 


0.040 


Elliptical Boundary 


12.15 


30 


0.998 



25 With respect to the interaction threshold boundaries considered, the results 

presented in Table 2 suggest that the elliptical boundary is the most appropriate to use in 
the interaction threshold model. Quasi-likelihood ratio tests given in (3.2.1) and (3.2.2) 
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are used to test the hypotheses of overall additivity and the hypothesis of the existence of 
the interaction threshold boundary in the experimental region, respectively. These results 
and the parameter estimates associated with the selected interaction threshold model are 
shown in Table 3. 

5 

Table 3 Parameter Estimates and Quasi-Likelihood Ratio Tests (Elliptical Boundary) 



Parameter 


Parameter 
Estimate 


Hypothesis Tested 


Test 
Statistic 


df 


p-value 


Po 


-21.85 


Additivity 


7.51 


6, 30 


<0.001 


Pi 


22.95 


Presence of Interaction 
Threshold Boundary 


6.78 


5, 30 


<0.001 


P 2 


27.05 










Pl2 


297.09 S 










A 


1.20 










B 


2.06 I 










C 


-1.74 










D 


-2.86 










E 


2.07 











* The value for f is 0.405, and the parameter estimates are associated with scaled 



10 doses. 

The hypothesis of additivity is rejected and the hypothesis that the interaction 
threshold does not exist is rejected, with p-values less than 0.001 for each test. It is 
interesting to consider the contours of the fitted dose response surfaces for this interaction 

15 threshold model. Since the isobologram presented by Gessner and Cabana (Figure 23) is 
associated with a 50% level of response, a comparison of the ED 5 o contours estimated 
from the interaction threshold model with the Gessner and Cabana isobologram is of 
interest, especially since the data set used in our interaction threshold model is an 
independently observed and much smaller data set. The estimated interaction threshold 

20 boundary plotted with the number of animals responding (out of 6) at each design point is 
graphically represented in Figure 25 and the contours associated with various EDioop's 
are found in Figure 26. 

The data set reported by Carter, et.al. (1988) that was used in the example is much 
smaller than Gessner and Cabana's experiment, and it includes many observations at dose 
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combinations where all or none of the animals responded, as illustrated in Figure 25. The 
similarity between Gessner and Cabana's ED 50 contour and that estimated by our 
interaction threshold model lends support to the concept of interaction thresholds. 
5. Discussion 

5 Gessner and Cabana's (1970) data suggests the existence of an interaction 

threshold for the combination of ethanol and chloral hydrate. The model presented here 
permits the estimation of the location of the interaction threshold boundary, and also 
supports such a finding. For more than two chemicals, the interaction threshold 
boundaries become multiple dimensional planes or surfaces, with increasing numbers of 

10 parameters needed to describe them. As the number of chemicals in the mixture 
increases, it becomes less likely that an investigator will be able to perform an 
experiment large enough to support the estimation of these parameters as well as those 
associated with the dose-response relationship. For mixtures with many components, it is 
likely that alternative approaches will need to be considered to account for these 

1 5 complexities with respect to interaction thresholds. 

The elliptical interaction threshold boundary raises some interesting points 
regarding interaction. Here, we rescaled doses for the initial purpose of imposing a 
continuity constraint in such a way that a point further from the boundary will have 
greater contribution to the interaction than a point that is near the boundary. There are 

20 many ways to accomplish such an effect, and the user has the flexibility to define the 
interaction effect thought to be most appropriate. As in the example presented, goodness- 
of-fit tests can be helpful in determining the appropriateness of the various choices 
considered. 

For studies of single agents, much has been written on the existence of dose 
25 thresholds (Lutz, 2001; Rroes, et.al, 2000; Sofuni, et.al 2000) and the development of 
models to accommodate the estimation of the threshold value. Gennings, et.al. (1997) 
developed additivity threshold models for drug/chemical combinations. We have added 
to the work on thresholds by introducing the concept of an interaction threshold. From 
our review of the statistical literature, the concept of an interaction threshold has not been 
30 considered. However, more recent studies on complex mixtures (El-Masri, Yushak, and 
Mumtaz 2002; Dobrev, Andersen, and Yang 2001) have also examined the assessment of 
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an interaction threshold in chemical mixtures through the use of physiologically based 
pharmacokinetic (PBPK) modeling. 
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